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Abstract: We apply the MC@NL0 approach to the process of heavy flavour 
hadropro duct ion. MC@NL0 is a method for matching next-to-leading order (NLO) 
QCD calculations and parton shower Monte Carlo (MC) simulations, with the fol- 
^ ■ lowing features: fully exclusive events are generated, with hadronisation according 

to the MC model; total rates are accurate to NLO; NLO results for distributions 
are recovered upon expansion in a s ; hard emissions are treated as in NLO compu- 
tations while soft /collinear emissions are handled by the MC simulation, with the 
same logarithmic accuracy as the MC; matching between the hard and soft regions 
is smooth, and no intermediate integration steps are necessary. The method was ap- 
plied previously to the hadroproduction of gauge boson pairs, which at NLO involves 
only initial-state QCD radiation and a unique colour structure. In heavy flavour pro- 
duction, it is necessary to include contributions from final-state QCD radiation and 
different colour flows. We present illustrative results on top and bottom production 
at the Tevatron and LHC. 
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1. Introduction 

The process of heavy flavour production in hadron collisions is a valuable testing- 
ground for perturbative QCD, since the high scale set by the quark mass should 
ensure that perturbative calculations are reliable. The prediction of cross sections 
and final-state distributions in heavy flavour production is also important for the 
design of collider experiments and new particle searches, since this process gives rise 
to irreducible backgrounds to many types of new physics. 

Up to now, the theoretical emphasis has been on next-to-leading order (NLO) 
calculations of total rates, single-inclusive distributions, and heavy quark-antiquark 
correlations, sometimes with resummation of higher-order contributions that are en- 
hanced in certain kinematic regions. However, for many purposes, such as studies 
of backgrounds to new physics, one needs a more complete characterization of the 
final state. This is provided by a Monte Carlo event generator program, which com- 
bines a calculation of the hard production process with a parton shower simulation 
and a hadronization model, to yield an approximate but realistic hadron-level event 
structure. 
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The problem with existing Monte Carlo event generators is that they are based on 
a leading-order (LO) calculation of the production process combined with a leading- 
logarithmic (LL) treatment of higher orders via the parton shower approximation. 
It has proved highly non-trivial to incorporate the benefits of NLO calculations into 
event generators, since the parton showers include parts of the NLO corrections, 
which should not be double-counted. On the other hand, the parton showers cannot 
be omitted, since they provide a reliable description of how final state hard partons 
evolve into QCD jets. Furthermore, any viable hadronization model operates on the 
multiparton states that are created by showering. 

A recent proposal for combining NLO calculations and parton showers is the 
so-called MC@NLO approach, introduced in ref. [1] (hereafter referred to as I). It 
is based on the highly successful subtraction method for NLO calculations. The 
basic idea is to modify the subtraction to take into account the terms that are 
generated by the parton shower. This results in a set of weighted LO and NLO 
parton configurations that can be fed into the parton showering generator without 
fear of double counting. Each weight distribution is well-behaved in the sense that it 
has no divergences or pathological tails that would lead to Monte Carlo inefficiency 
However, in order to reproduce the NLO corrections fully, some of the configurations 
have negative weights. Event unweighting can still be achieved efficiently, if desired, 
by generating a small fraction of 'counter-events' that contribute with equal but 
opposite weight to events in all distributions. 

The MC@NLO method was worked out in detail in I for processes in which, at 
the Born level, there are no coloured partons in the final state. An important example 
is gauge boson pair production, for which a wide range of MC@NLO predictions were 
presented there. To deal with the process of heavy quark production we must take 
into account QCD radiation from final-state partons, in this case the heavy quarks 
themselves. A further new complication is the possibility of different colour flows. 

We shall see that no difficulties of principle arise from these complications. The 
only new task, albeit a laborious one, is to calculate precisely what the shower Monte 
Carlo is doing at NLO, in order to compute the modified subtractions correctly. In 
our case we use the HERWIG shower Monte Carlo [2] , the relevant features of which 
are summarized in Appendix A. 

In the following section we review the main features of the MC@NLO approach. 
Then, in sect. 3, we discuss the partonic processes that contribute to heavy flavour 
production in standard Monte Carlos and in MC@NLO, and we define the kinematic 
variables for the corresponding 2^2 and 2^3 processes. In order to compute 
the Monte Carlo subtraction terms that form the basis of the MC@NLO method, 
we have first to relate these variables to those used in the HERWIG program. This 
is done in sect. 4. Next we write down, in sect. 5, the approximate 2^3 particle 
production cross sections generated by the HERWIG parton showering algorithm, 
which are then used to construct the Monte Carlo subtraction terms for heavy flavour 
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production. Technical details of this procedure are given in Appendix B. Inserting 
the subtraction terms in the formulae reviewed in sect. 2 enables us to generate 
parton configurations that can be fed into the shower Monte Carlo without any 
double counting of NLO contributions. In order for the HERWIG Monte Carlo to 
operate correctly we have to assign a colour flow to each configuration; this and other 
details of the implementation are explained in sect. 6. 

We present the predictions of MC@NLO for top quark production in sect. 7. 
The case of bottom is much more involved than that of top. Problems affecting b- 
physics simulations with standard Monte Carlos are reported in sect. 8.1; in sect. 8.2, 
we discuss some of the features of MC@NLO in b production, and in particular the 
treatment of large logarithms of p T /iri, the ratio of the quark transverse momentum 
to its mass; in sects. 8.3 and 8.4 we present MC@NLO predictions for bb correlations 
and single-inclusive distributions respectively, and compare them to HERWIG and 
NLO results. Finally, conclusions and future prospects are presented in sect. 9. 



2. Review of MC@NLO approach 



The MC@NLO formalism is defined in eq. (4.22) of I, which we denote by eq. (1.4.22). 
We rewrite that equation in the following, fully equivalent, form: 
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The quantities da/dO, I MC (0,3), and I h 



d<p2dx 

(O, 2) appearing in eq. (1.4.22) have been 



replaced here by .Fmc@nlo (the MC@NLO generating functional), (the MC gen- 
erating functional when starting from a 2 — > 3 hard subprocess), and J^^c (the MC 
generating functional when starting from a 2 — > 2 hard subprocess) respectively. 
This renders more transparent the fact that the dependence upon the observable 
O in eq. (1.4.22) is only formal, and that the MC@NLO generates events without 
reference to any observable. We refer the reader to I for the definitions of all the 
terms appearing in eq. (2.1); however, the precise details will not be relevant here. 
In what follows, we shall limit ourselves to describing the basic features of eq. (2.1), 
and their role in the implementation of heavy flavour production in MC@NLO. 

We start by recalling that the definition of MC@NLO is derived from the ex- 
pectation value (O), computed at the NLO, of a generic observable O, eq. (1.4.19). 
This can be read off from eq. (2.1) simply by removing the MC subtraction terms 
c?E a6 | MC , and by replacing with O(n), where O(n) is the observable O computed 
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in an n-body final-state configuration (n = 2,3). The NLO expression for (O) is an 
integral over the momentum fractions x\ and X2 of the incoming partons, and three- 
body phase-space variables 3 . Each point (xi,x 2 ,0 3 ) in the integration range thus 
corresponds to a 2 — > 3 kinematic configuration (called EI hereafter). Furthermore, 
a definite 2-^2 configuration is chosen (called § hereafter), according to a mapping 
which we denote by Va^s- The weight associated with EI is given by the real-emission 
matrix element; the sum of all the remaining contributions to the NLO cross section 
(namely, the Born term, the virtual term, the soft and collinear counterterms, and 
the finite remainders of the initial-state collinear singularity subtractions) constitutes 
the weight associated with §. One can prove that it is always possible to cast any 
expectation value (O) in this form (see I, sects. 4.4 and A. 4), through a formal pro- 
cedure that we call event projection. It should be clear that event projection does 
not imply any approximation, and that its specific form depends on the choice of 
Vu^n- In the context of a pure NLO computation, this choice is arbitrary, and its 
freedom has been used in the past [3] to improve the convergence of the numerical 
integration procedure. On the other hand, when defining an MC@NLO it is the MC 
itself that dictates the form of Ph-*§- 

We suppose now to have chosen a map Vu-,s, and to have performed event 
projection on an NLO cross section. For each point (xi,X2,(f>3) we get a pair of 
kinematic configurations, EI and S. Instead of using these configurations for defining 
the observable O, as in the NLO computation of (O), we feed them into an MC, 
where they are treated as initial conditions for shower evolution. This corresponds 
to defining the MC@NLO as in eq. (2.1), except for the MC subtraction terms cE ab | MC 
which are omitted. However, this naive attempt fails (see I, sect. 3.3.1). Basically, 
when evolving § configurations the shower reproduces some of the EI configurations, 
which are therefore double counted. The idea of modified subtraction, upon which 
the MC@NLO approach is based, is to subtract these double-counted configurations 
at the level of short-distance cross sections. This is the role of the MC subtraction 
terms d£ a6 | MC which appear in eq. (2.1); they come in pairs, since they have to 
account for both the emission and the non-branching probabilities in the MC. The 
MC subtraction terms act as local counterterms in eq. (2.1), and this implies that 
the weight distributions for H and § configurations (the terms multiplying ^fc an d 
^ic respectively, see also eq. (1.4.23) and eq. (1.4.24)) are separately convergent, 
thus allowing event unweighting as is customary in MC simulations. 

The MC subtraction terms are obtained by formally expanding the MC results to 
the first non-trivial order in a s , which corresponds to an EI configuration. Typically, 
d^ablmc has the form of a hard, 2^2 cross section, times a kernel which describes 
(the first) parton branching. The shower algorithm fully specifies how to determine 
the EI configuration, given the hard § configuration and the values of the showering 
variables. Thus, it implicitly defines a map between the EI and § configurations. We 
choose Vw^s, used in event projection, to coincide with this map. Notice that this 
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is possible only because the shower algorithm is independent of the hard process. 

In summary, the main steps that have to be taken in order to construct an 
MC@NLO are the following. 

i) Determine Vw^s- 

ii) Write the NLO cross section for the relevant production process, and perform 
event projection on it, using the map Vm^n found in %). 

Hi) Define the MC subtraction terms, and insert them into the expression for the 
NLO cross section. 

Of the three steps above, only ii) depends on the process in a non-trivial way. The 
implementation of event projection requires a detailed knowledge of the formalism 
adopted to write the cross section to the NLO accuracy; an explicit example of such 
a procedure has been given in I, but different prescriptions are clearly possible. On 
the other hand, step i) is strictly process- independent, and thus Vu^s can be deter- 
mined once and for all. However, Vw^s depends on the particular shower algorithm 
adopted, and therefore different MC's define different Vw^s maps. Finally, step Hi) 
is process dependent, but only through the hard 2^2 cross sections that appear in 
a factorized form in the MC subtraction terms, i.e. at the LO level, which is fairly 
simple to deal with. The part of the MC subtraction term which describes the first 
branching depends only on the shower algorithm, and can therefore be studied with 
full generality. 

The definition of the MC subtraction terms is in fact one of the main goals 
of the present paper; final-state emissions are considered here for the first time, 
since they were not relevant to the production of vector boson pairs considered in 
I. Also, the formulae presented in I for initial-state emissions will be generalized 
here to account for colour structures more complicated that those of I. These results 
will allow an almost straightforward definition of the MC subtraction terms for any 
production process. A subtle point concerns the interplay between initial- and final- 
state emissions, and its impact on the definition of S-event contribution to MC@NLO. 
This issue is discussed in Appendix B. 

3. Heavy quark production 

3.1 Contributing processes 

In the MC@NLO approach, the Monte Carlo is not involved in the generation of the 
hard process which, apart from the modified subtraction, is treated as in standard 
NLO codes. It follows that the partonic production processes that we need consider 
are 

qq^QQ, gg -»• QQ (3.1) 
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at 0{a: 



:|) and O(al), and 



qq -> ^ -> #g -> QQg, gq -> QQg 



(3.2) 



at 0(«g). As discussed previously, the C?(a|) contributions in eq. (3.1) generate 
(some of) the configurations in eq. (3.2) through parton showering, but MC@NLO 
is defined in such a way to avoid any double counting. The processes in eq. (3.1) are 
traditionally classified as flavour creation (FCR hereafter) to distinguish them from 
the so-called flavour excitation (FEX hereafter) processes 



where charge-conjugate processes are understood to be included. In the case of FEX 
processes, it is implicitly assumed that the relevant heavy flavour is already present 
in the parton distribution functions (PDFs) of the incoming hadrons. 

In standard Monte Carlo, both FCR and FEX matrix elements are used in the 
generation of the hard processes. 1 This is not the case in NLO computations, since 
the partons which initiate the hard processes are treated as massless, understanding 
that their actual mass is negligible with respect to the hard scale of the process, in 
this case the heavy quark mass m. This forbids the presence of heavy quarks in 
the initial state. Furthermore, if we use NLO matrix elements in the generation of 
the hard process, as in MC@NLO, the distinction between FCR and FEX becomes 
ambiguous. For example, the NLO FCR process gg — > QQg has a contribution 
from initial-state quasi-collinear gluon splitting, g — > QQ, followed by the LO FEX 
process Qg — > Qg. Since initial-state quasi-collinear gluon splitting forms part of 
the evolution of the PDFs of the incoming hadrons implemented through parton 
showers, we are in danger of double-counting if we include both LO FEX and NLO 
FCR processes in MC@NLO. On the other hand, this example explains why FEX 
processes are considered in standard Monte Carlo: they allow one to include some 
of the features that could not be included by shower evolution initiated by FCR 
processes. We also point out that the process gg — > QQg, which we include in 
the MC@NLO, generates certain kinematic configurations that can be equivalently 
produced starting from a hard gg — > gg process followed by final-state gluon splitting 
g — > QQ. The gluon splitting mechanism (GSP hereafter) is known to be important 
in bottom production, giving the dominant contribution to those configurations in 
which b and b are close in phase space (which is important not only for the study 
of correlations, but also in 6-jet physics), and non- negligible contributions to single- 
inclusive distributions in the intermediate-^ region. 

Flavour excitation and gluon splitting pose a couple of non-trivial problems in 
standard Monte Carlo event generators. 2 In FEX, the presence of a heavy quark in 
X FEX is not relevant to top production at realistic collider energies. 



2 For a discussion of b physics as described by Monte Carlo parton shower programs, see e.g. 
refs. [4, 5]. 



qQ -> qQ, qQ -> qQ, gQ -> gQ 



(3.3) 
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the initial state implies a direct dependence upon its PDF; in the small-p T region, 
i.e. close to threshold, it is unlikely that the treatment of the backward evolution 
performed by the parton shower will be compatible with what was implemented in 
the PDF evolution. Furthermore, the t-channel singularity in the matrix elements re- 
quires the implementation of a cutoff (possibly an effective one, see app. A. 1.2) which 
prevents the generation of p T = events. Both features result in predictions which 
have a certain degree of implementation dependence. As far as GSP is concerned, 
the corresponding events are usually obtained by considering pure QCD events (i.e., 
with no heavy quark involved in the hard process, the generation of which also re- 
quires a p T cutoff), and selecting those in which at least one QQ pair is obtained 
through showering, a very inefficient procedure from the statistical point of view. 
More details will be given in sect. 8.1. 

Worst of all, both FEX and GSP processes are well defined only in the case 
of large transverse momenta of the heavy quark. Their extrapolation to the low 
transverse momentum region can at best be considered a very rough model of higher- 
order heavy flavour production processes. 

The implementation of heavy flavour production in MC@NLO helps to avoid 
some of these problematic features of the FEX and GSP processes. In the processes 
considered in NLO computations, eqs. (3.1) and (3.2), there is no such thing as 
direct dependence upon the heavy quark PDF, and singularities are cancelled at 
the level of short-distance cross sections, without the need to introduce unphysical 
parameters acting as cutoffs. Furthermore, one generates gluon splittings with 100% 
efficiency, since they are included in the matrix elements. Clearly, the presence of 
negative weights in MC@NLO degrades the statistics, but (at least with the fraction 
of negative weights we find in our implementation) the overall efficiency is still larger 
than that required to obtain a useful gluon-splitting sample with the usual Monte 
Carlo method. 

However, an MC@NLO based on the hard processes of eqs. (3.1) and (3.2) cannot 
account for all the contributions generated by FEX and GSP. For example, MC@NLO 
does not include the diagrams in which a gluon, rather than a heavy quark, is emitted 
in the first backward branching of the heavy quark entering a FEX hard process. 
Similarly, MC@NLO does not include those GSP diagrams in which the gluon emits 
another gluon before splitting. Examples of omitted contributions are shown in fig. 1. 

The absence of such contributions is not surprising: the terms missed by the 
MC@NLO are part of those relevant to the region of large heavy-quark transverse 
momenta (large p T ), where the quark mass m no longer sets the hard scale. In 
this region, the NLO computation is not expected to give sensible results, since 
large logarithms log(p T /m) need to be resummed to all orders. For single-inclusive 
observables, techniques are known for resumming these large logs (to NLL accuracy), 
and for matching the resummed result to the NLO one (for example, FONLL [6]). 
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Figure 1: Graphs that give rise to enhanced higher-order terms, and are not included in 
the MC@NLO implementation. Shaded areas represent almost collinear branchings due to 
parton showers. 

However, there is at present no known way of extending such a procedure to the fully 
exclusive case that would be needed in order to include it in MC@NLO. This problem 
is relevant, for example, to bottom production at the Tevatron. In sect. 8.2 we shall 
show that MC@NLO gives in practice a sensible answer also in this case, since these 
enhanced high-p T effects are not so large in the kinematic region of practical interest. 

3.2 2 — > 2 processes 

In this and in the following section, we introduce the notation that we shall use in 
order to describe the hard processes in eqs. (3.1) and (3.2). 

We denote variables relating to the 2 — ^ 2 FCR processes of eq. (3.1) by barred 
symbols: pi,p2 represent the incoming (massless) momenta and k\,k2 the outgoing 
heavy quark and antiquark (mass m) momenta. The momentum fractions of the 
incoming partons are denoted by x±,2, i.e. 

pi = x x P\ , p2 = x 2 P 2 , (3.4) 

where are the beam momenta, and the invariants by 

s = 2pi ■ p 2 , t — -2pi ■ ki , u = -2pi ■ k 2 . (3.5) 

Then s + i + u = and s = X1X2S, where S is the overall cm. energy squared. In 
the cm. frame of the 2 — » 2 process, we can write 

p lj2 = £(l,0,0,±l) , h,2 = (E,±k T ,0,±k L ) , (3.6) 

where 

s = AE 2 , t = -2E{E -h) , u= -2E(E + k L ) ■ (3.7) 
We also introduce the cm. scattering angle 9 and heavy quark velocity 

=y /l- Am 2 /s , (3.8) 

so that 

t = -is(l-/3cos£) , u = -|s(l + (3cos9) . (3.9) 
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The hard scattering Born cross section is then 



da = dx! dx 2 /i Hl) (x!)/ 6 (H2) (x 2 ) da ab , (3.10) 

ab 

where 

da ab = Mabd(p2 = -^MabdcosO, (3.11) 

Ai a b being the spin- and colour-averaged Born matrix element squared for the process 
ab — > QQ, times the flux factor: 

— 4 iV 2 -l 1 (\ iu m 2 \ 

-T-i 4 N 1 (u t 1 s 2 \ fl tu 2m 2 2m 4 \ , oio . 
3.3 2-^3 processes 

For the 2 — > 3 FCR processes of eq. (3.2) we use unbarred symbols: p±,P2 for the 
incoming (massless) momenta, k±, k 2 for the outgoing heavy quark and antiquark 
momenta, respectively, and k for the momentum of the emitted (massless) parton. 
We denote the momentum fractions of the incoming partons by x^2, i.e. 

pi = xiPi , p 2 = x 2 P 2 , (3.14) 

and the invariants by 

s = 2p l ■ p 2 , ti = -2pi ■ ki , t 2 = -2p 2 ■ k 2 , Mi = -2p l ■ k 2 , m 2 = -2p 2 • h , 

(3.15) 

so that s = X\X 2 S. We also introduce 

vi = —2pi ■ k , w 2 = -2p 2 • k , (3.16) 

and 

Wi = 2ki ■ k , w 2 = 2/c 2 • /c , (3-17) 
but these are not independent variables since 

s + ti + Ui + Vi = s + t 2 + u 2 + v 2 = , 

s + ti + u 2 — w 2 = s + t 2 + -Ux — w\ = . (3.18) 
This notation is summarized and compared with that of ref. [3] in table 1. 
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Label 


Invariant 


Ref. [3] 


Relation 


s 


2v-i ■ z>o 


s 






—2pi ■ k\ 






to 


— 2vo ■ ko 


HI 




^ 1 


— 2oi • ko 


yi 




Ll Z 


— 2t>o ■ ki 


do 




L 1 


—2pi ■ k 


T, h. 


— S — t\ — U\ 


v 2 


-2p 2 ■ k 


u k 


—s — t 2 — u 2 


Wi 


2ki ■ k 


Wi 


S + t 2 + Mi 


w 2 


2k 2 ■ k 


w 2 


S + ti + U 2 




(h + k 2 f 


s 2 


S + Vi + V 2 



Table 1: Notation for 2 — > 3 kinematics. 



4. Relating NLO and MC kinematics 

As discussed in sect. 2, the implementation of MC@NLO requires a detailed knowl- 
edge of the relation between 2^3 (H) and 2 — > 2 (§) kinematic configurations. 
Also, for the integration of the MC subtraction terms over the three-body phase 
space, as prescribed in eq. (2.1), we need to express the MC showering variables 
in terms of the phase-space variables. We solve both problems in the same way, 
presenting the 2-^2 invariants and showering variables as functions of the three- 
body invariants. The latter will eventually be computed by means of the three-body 
phase-space variables; however, in this section we shall not introduce any explicit 
phase-space parametrization, and the results are therefore fully general. 

A shower Monte Carlo program generates arbitrarily complicated multiparton 
configurations by an iterative process of parton emission, starting from the coloured 
external lines of a hard subprocess. In the hard subprocess these external lines 
are treated as being on mass-shell, but the showering process may drive them off 
the mass shell, and the kinematics have to be adjusted to restore energy- momentum 
conservation. This process is called momentum reshuffling; the way it is implemented 
in the HERWIG program is described in app. A. 2. 4. The relationship between 2 — > 
2 and 2^3 kinematics discussed in the following subsections depends upon the 
particular implementation of momentum reshuffling. 

4.1 Relating 2^2 and 2^3 kinematics 

The relationship between the variables of the 2 — > 2 hard process, and the 2-^3 
variables which result after one parton emission is not simple, owing to momentum 



11 



reshuffling. In particular, since initial- and final-state showers are treated differently 
in the reshuffling, it depends on whether the parton emission is from an incoming or 
outgoing parton. 

We consider the case of FCR with 2-^2 momenta pip 2 — > kik 2 and 2^3 
momenta p\p 2 — > k\k 2 k. The way in which we relate the kinematics of the 2-^2 
and 2 — > 3 processes is by solving for the 2^2 invariants s, i, u and momentum 
fractions X\, 2 in terms of the 2-^3 invariants s, £1,2,^1,2 an d momentum fractions 
Xi >2 . In terms of the cm. frame variables defined in eq. (3.6), this amounts to finding 
Xi )2 , E, and k h . 

It should be noted that the kinematics used in parton shower Monte Carlos 
generally involve cutoffs that operate like effective light quark and gluon masses. We 
can ignore these cutoffs in computing the Monte Carlo subtraction terms since they 
only give rise to power-suppressed corrections, comparable to hadronization effects, 
in physical distributions. This point was discussed in I, sect. B.3. 

4.1.1 Final-state emission 

Suppose a gluon is emitted from the heavy quark. The final state is then formed 
by the heavy quark jet (i.e. the heavy quark plus the emitted gluon) and the heavy 
antiquark. The three-momenta of the antiquark and of the heavy quark jet are 
rescaled by a common factor in order to restore energy conservation, according to 
the prescription described in appendix A. 2. 4. The heavy antiquark has momentum 
(in the hard process cm. frame) 

k 2 = (V m 2 + a 2 k 2 , —ak T , 0, — ak h ) (4.1) 

where k 2 = k 2 + k 2 = E 2 — m 2 and a is the rescaling factor, given by the energy 
conservation constraint, 

2E=^m 2 + (ak) 2 + \/{k + h) 2 + (ok) 2 . (4.2) 
Hence for final- state emission from the heavy quark Q we find 

E = Wi, A; L = ^f^i)A, (4.3) 
V s - wi J f3 2 

where (5 2 is the velocity of the heavy antiquark in the heavy quark-antiquark cm. 
frame, 

= y / l-Asm 2 /(s-w 1 ) 2 . (4.4) 

For emission from the heavy antiquark, the labels 1 and 2 are interchanged, so 
the relation between the 2^2 and 2^3 invariants depends on which outgoing 
parton emits the gluon. We label the 2^2 invariants as i Q) Iq, etc., according to 
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which parton is the emitter. Then we have explicitly 



S Q = 


S Q = 


S 


h = 


-is 


1 - 


h = 


-is 


1 - 




—s — 





h — Ul\ P 



s -w 1 J p 2 
h — u 2 \ (3 



where 



s - w 2 J Pi 
uq = -s-i 



Pi = a/1 - 4sm 2 /(s - w 2 )' 
0, £1,2 -> t, ui :2 



u 



we have tq q 



(4.5) 

(4.6) 
t and 



Note that in the soft limit (wi j2 
Uq q — > u as expected. 

The incoming parton momenta are not affected, so pi = pi, p 2 = p 2 and Sq q = s 
Thus the incoming momentum fractions are unchanged by final-state emission: 



Xif = Xi , x 2f =x 2 



(4.7) 



The formulae presented here can also be used for the branchings of massless 
partons; one simply lets j3 — > 1, Pi — > 1, p 2 — > 1. In such cases, collinear limits must 
also be considered. When a gluon is emitted collinearly by the heavy quark (t 2 — > i, 
Ui — > u, ti + Vi — > t, u 2 +v 2 ^ u) we have Iq — > i, uq — > -u. For a collinear emission 
by the heavy antiquark (ti — > t, m 2 — > m, ^2 + ^2-^ ?, «i + ^i — > m) we have tg — > t, 

UQ -»■ «. 

4.1.2 Initial-state emission 

In the case of emission from the incoming partons, the invariant mass of the heavy 
quark pair is not changed by momentum reshuffling, so 



s — (ki + k 2 ) 2 = s + Vi + v 2 . 



(4.8) 



However, the pair receives longitudinal and transverse boosts, as described in app. A. 2. 4. 
To relate the other invariants, we consider the heavy quark momentum difference. 
In the 2 — > 2 cm. frame, it is 



fci - fc 2 = 2(0, fc T , A) ■ 



(4.9) 



Since this has no energy component, the effect of the longitudinal boost is simply to 
rescale the longitudinal component by the boost factor 



1L 



1 + 



(h + k 2 )\ 
(h + k 2 f 



(4.10) 
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The transverse boost does not change the longitudinal component, so 



(h - k 2 ) L = lL (h - k 2 ) L = 2 7z A . (4.11) 
We can extract the longitudinal components using the combination 

■ r - 1, '\ - XU '- 1 m.n.n.n. (4.12) 

y/XiX 2 S 

Hence 

ib L = _(^i-^)-( fc i-^) > (4.13) 

where 



/-, [{X2P1 -xm) • (fci + k 2 )] 2 

7l = \1-\ 71 — — pva • ^ 4 - 14 ^ 

xix 2 s(A;i + ^2) 

Expressing everything in terms of invariants, we therefore have for initial-state emis- 
sion 

J7 1 / , , T ^ X 2 (h - Ux) + X!(t 2 - u 2 ) 

E = ±Vs + v 1 +v 2 , k^ = E— — — , (4.15) 

2S - £i£ 2 V \V 2 jS l 

where for future reference we have defined 

1 fs + v 2 s + v x \ 

x± = - xi ± x 2 . (4.16) 

2 \ s s J 

Note that in this case the result is independent of which incoming parton emits. 
However, for consistency we label t = t + or t- according to whether parton 1 or 2 is 
emitting. Thus we have explicitly 



s± = s + v 1 + v 2 
t± = -|(s + Vi+ v 2 ) 



_ x 2 {t 1 - Ui) + Xi(t 2 - u 2 ) 
2s^Jx\ - xix 2 viv 2 /s 2 
u± = -s± - i± . (4.17) 



Again in the soft limit (t>i,2 — > 0, ti j2 — > t, u\^ 2 — > u) we have t± — > t and w± — > w 
as expected. In the case of collinear emission from leg 1 (t>i — > 0, £2 - ► u 2 — > w, 
ii + tui — > t, Mi + cj 2 — > m) we have i± — > £, ti± — > w. For a collinear emission from leg 
2 (f ! — > 0, £1 — > t, Mi —> u, t 2 + uj 2 ^> i, u 2 + c<Ji — > m) we also have t± — > £, w± — > -u. 

To relate the incoming momentum fractions, we note that in this case eq. (4.8) 
applies, and therefore instead of eq. (4.7) we have 

- - s + vi+v 2 

xux 2i = x ± x 2 . (4.18) 

s 

The values of x Vl and x 2i separately depend on the momentum reshuffling scheme, 
as explained in app. A. 2. 4: 
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• In the p-scheme the longitudinal momentum of the heavy-quark pair is pre- 
served, which means that 

S + V2 S + Vl I A 1(V» 

xii-x 2i = xi x 2 , (4.19) 

s s 

and hence 

x u = x_ + \jx\ - x\XiV\Vil s 2 , x 2i — xu - 2x- , (4.20) 

where x± are given by eq. (4.16). Note that eq. (4.20) coincides with eq. (I. A. 47); 
this is to be expected, since the details of the hard process are irrelevant in the 
determination of the x's. 

• In the y-scheme the rapidity of the heavy-quark pair is preserved, which implies 

Xu x 1 (s + v 2 ) 



X2i X 2 (S + Vi) ' 

and hence 



(4.21) 



(s + V! + v 2 )(s + v 2 ) _ (s + v 1 + v 2 )(s + v 1 ) 

x li = x 1 * — r , X 2i = X 2 \ — r , 4.22 

y s(s + vi) y s(s + v 2 ) 

which corresponds to eq. (I. A. 42). 
4.2 Relating HERWIG variables to invariants 

When relating the HERWIG showering variables to the kinematics of the 2^3 hard 
process, we must again take careful account of which quantities are preserved under 
momentum reshuffling. 

4.2.1 Final-state emission 

For emission of a gluon from the outgoing heavy quark in the FCR processes, the 
HERWIG variables are the angular variable £ = k ■ ki/k^k^ and the energy fraction 
z = k®/Eo = 1 — k /Eo, all energies being evaluated in the showering frame where 
Eq = —t/2 or —u/2 as discussed in app. A. 2.1. The invariant quantities are the jet 
virtuality 

(Jfei + k) 2 - m 2 = Wl = 2z(l - z)£E% , (4.23) 
and the "+" momentum fraction of the gluon with respect to the jet axis, 

Cl = 7r A^ = (1 _, ) I±MM, (4.24) 
(h + k) ■ n 2 i+ft 

where (3\ is the heavy-quark jet velocity in the showering frame 



fc = \ l-( Wl +m 2 )/E 2 , (4.25) 
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and n 2 is a lightlike vector along the direction of the heavy antiquark in the heavy 
quark- ant iquark cm. frame: 



n 2 = h- ^"2^(1 - &)(pi +P2) , (4.26) 

where /3 2 is given by eq. (4.4). Inserting eq. (4.26) into eq. (4.24), we find 

c = (s + wi)w 2 + (s- wi)[(wi + w 2 )(3 2 - Wi] 

^ (s- Wl )/3 2 [( S + Wl ) + ( S - Wl )/? 2 ] • 1 • } 

Solving eqs. (4.23), (4.24) and (4.27) for z and £ with £g = -t Q /2 as given in 
eq. (4.5), we obtain z*q and corresponding to emission from the heavy quark 
with the t-flow colour structure as defined in app. A. 1.3. Similarly, solving with 
Eq = —uq/2, as appropriate to the -u-flow, gives Zq^ and £q\ We find that 

z§ = l- ACi ^— , (4.28) 

where Iq = Iq, uq for I — t,u. Interchanging the labels 1 and 2 and using Iq instead 
of Iq gives Zq' w ^ and £^'"' ) , corresponding to emission from the heavy antiquark. 

Note that the transverse momentum of the emitted gluon relative to the jet axis 
is k T where 

^ = [(l-Ci)^i-G™ 2 ] , (4.30) 

and thus k T ~ \/2z(l — z)Q at small values of £ and m 2 /Eq, where Q = E ^ is the 
HERWIG evolution variable. 

4.2.2 Initial-state emission 

In the case of an initial-state jet, the jet axis coincides with the beam axis and so 
the jet invariants are simpler than in the final-state case. For emission from parton 
1, the jet virtuality is 

(p 1 -k) 2 = v 1 = -2±^-Zl%, (4.31) 
and the "+" component of the gluon momentum is 3 

k ■ p 2 v 2 



Pl'P2 



Ul-z)(2-0- (4.32) 



Solving eqs. (4.31) and (4.32) for z and £ with E 2 = —t+/2 or — u + /2 as given in 
eq. (4.17), we obtain z+' u ^ and corresponding to emission from incoming parton 

3 Note that eqs. (4.31) and (4.32) are equivalent to eq. (I. A. 67) and eq. (I. A. 68), if we set 
E 2 Q =M^ w /2. 
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1 with these initial conditions. In the case of emission from an initial-state gluon, 
we also have the possibility that E 2 = s + /2 (see app. A. 1.3), and we denote the 
corresponding variables by and Then we can write the solutions explicitly 

as 



£0 



ft 



Vl 



Vl 



1 + 



V2 



1 + 



V2 



s(l-z. 



(4.33) 
(4.34) 



where l + = s+,i+, u + for I — s,t, u. 

For emission from incoming parton 2, the variables V\ and v 2 are interchanged. 
We denote the corresponding solutions by and ^ . 

4.3 Dead regions 

The parton shower initial conditions imply that gluon radiation is confined to an- 
gular regions (cones) specified by the colour flow. This represents an approximate 
treatment of the destructive interference due to colour coherence. The dead regions 
outside the cones can be found from the above mapping of the HERWIG shower 
variables onto invariants. 

4.3.1 Final-state emission 

For final-state jets the kinematic region available for gluon emission is 

m 



l-y/l-m*/E*<Z<l, 



< z < 1 



(4.35) 



However, the radiation in the forward direction is suppressed dynamically for £ < 
m 2 /Eq and therefore in HERWIG the region used (neglecting the gluon effective mass 
cutoff) is 



m 2 /El < £ < 1 , 



m 



< z < 1 . 



(4.36) 



The scale E is given by E 2 — — t/2 — s(l — (3 cos 9) /A or — u/2 or s/2, depending on 
the colour flow (see app. A. 2.1). Thus all cases can be represented by the situation 
for Eq = —t/2, with cos^^ — cos 9 when Eq = —u/2 and cos# = —1/(3 ~ —1 when 
E 2 = s/2. 

It is convenient to express the jet regions in terms of the Dalitz plot variables of 
the QQg final state, 



x Q = 2k x ■ (pi +p 2 )/s = -(ti + u 2 )/s = 1 - w 2 /s 
xq = 2k 2 ■ (pi +P2)/s = -(ui + t 2 )/s = 1 - wi/s 

X g = 2k ■ (pi +P 2 )/S = -(Vl + V 2 )/S = 2 - Xq - Xq . 



(4.37) 
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Then for emission from the heavy quark we have 

wi = s(l -Xq) , 

Ci 



x a + ( x 1 + X Q ~ x q)/(^ x Q^) 



2-x q (1-(3 2 



(4.38) 
(4.39) 



where W\ and Ci are related to the HERWIG variables £ and z by eqs. (4.23) and (4.24). 
For emission from the heavy antiquark, the variables xq and xq are interchanged. 
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Figure 2: Dalitz plot and jet regions for final-state emission when s = M§, m = 5 GeV, 
cos 9 =0.5 (dotted), (dot-dash), -0.5 (dashed) and -1 (solid). 



The boundaries of the quark and antiquark jet regions of the QQg Dalitz plot 
are shown in fig. 2. Note that there is an overlap in the soft region for cos^ < 0. The 
boundary of the physical region is also shown (short-dashed). We see that there is 
a dead region in which hard, non-collinear gluon emission is missed by the HERWIG 
shower algorithm, and also narrow near-collinear dead regions, where emission is 
forbidden in order to simulate the dynamical suppression of collinear emission from 
the heavy quarks, as discussed above. 

4.3.2 Initial-state emission 

In the case of an initial-state jet, the allowed region is £ < z 2 since now z = E /p® 
[7]. The conventional variables are x, y which give the heavy diquark mass Mqq and 
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the emission angle 9* of the gluon in the partonic cm. frame. For emission from 
incoming parton 1 we find: 

^M| s / S =l + HL±^, vsm<r (4.40) 

W S y Vi + v 2 

For emission from parton 2, v i and v 2 are interchanged, i.e. y — > — y. Inserting £ = 
in eqs. (4.31)-(4.32) then gives the boundaries of the jet regions in the x, y plane. 




0.2 0.4 0.6 0.8 1 

X 



Figure 3: Phase space and jet regions for initial-state emission when s = M§, m = 5 
GeV, cos£ =0.5 (dotted), (dot-dash), -0.5 (dashed) and -1 (solid). 

The jet boundaries are shown in fig. 3. Note that there is again an overlap in the 
soft region for cos(9 < 0. There are no collinear dead regions, because the incoming 
partons are treated as massless. 

5. MC cross sections expanded to NLO 

In this section, we present the cross section that we would obtain by keeping the 
first non-trivial order in the a s expansion of the HERWIG Monte Carlo result. This 
quantity, which we denote by da\ MC , is directly related to the MC subtraction terms 
that enter the definition of MC@NLO. We shall not give here the rather technical 
details of the construction of the MC subtraction terms, which we report in app. B. 
We shall limit ourselves to highlighting the main differences with respect to I, which 
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da 



result from the more complicated singularity and colour structure of the heavy flavour 
cross section. 

The result we are seeking can be written as follows: 

= EEE^L. (») 

ab L I 

where the first sum in eq. (5.1) runs over parton processes. The index L runs over 
the emitting legs and, consistently with sect. 4.1.1 and 4.1.2, it assumes the values 
+ , — , Q, and Q. The index I runs over the colour structures, and it can take the 
values s, t, and u (see sect. 4.2.1 and 4.2.2). Using the same formal expansion as in 
eq. (I. A. 58), we obtain 

1 



da 



(+,0 

ab 



da 



(-1) 

ab 



fi H ^ulz^\x 2i )da^ 
y^H)fb iH2 \^/z^)da^ 



(0 



dX\i d%2i i 



dX\i d%2i i 



da 



da 



(Q,i) 

ab 

(Q,i) 

ab 



e^f)fl H2 \-2 f )da^ 



dx\ f dx 



MC 



J a 



(x 2f ) d<x 



ab 



MC 



if ax 2 f 



dx±f dx 2 f 



(5.2) 

(5.3) 

(5.4) 
(5.5) 



In the case oil — s, the sum of eqs. (5.2) and (5.3) coincides with eq. (I. A. 58). The 
cases I — t,u, and eqs. (5.4) and (5.5) were not considered in I, since they correspond 
to colour flows not relevant to gauge boson pair production, and emissions from 
strongly-interacting final-state partons. 

The short-distance cross sections da in eqs. (5.2)-(5.5) have a form similar to 
eq. (I. A. 63) and eq. (I. A. 64), namely a factor depending on the HERWIG showering 
variables £ and z, times a Born cross section. As discussed above, £ and z, and the 
2^2 invariants s, i and u entering the Born cross sections, must be expressed in 
terms of the 2^3 integration variables used in the NLO code. In the case of emis- 
sions from both the initial- and final-state partons, we have four pairs of functional 
relations between £ and z and the NLO integration variables, corresponding to z^} , 

We have the following non-vanishing contributions: 
• qq initial state 



and z%' (and analogously for £), given by eqs. (4.33,4.34) and (4.28,4.29). 



da 



(+,*) 
qq 



2tt 



%dzfP^(zf)da q ^([zff-if) 



da^ t] 



da [Q -' l) 



da, 



(+,*) 



qq 



dt 



MC 
(*) 



,(*) 



2tt 

so 



■dzVp^d^e (i-$)e Uz$y 



da, 



(Q,t) 
qq 



da, 



(Q,t) 
qq 



(• 



Q 



(t) 



(*) 



4'0 
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• qg initial state 




where da g g are the colour t- and u-flow contributions defined in eq. (A. 6); 
• gg initial state 




Note that in each equation the Born cross section da a b or daj has to be computed 
using the relevant definitions of s, i and u; more details are given in app. B. The 
argument of a s used in the computation of these cross sections is the same one as 
that used in the computation of the NLO short- distance cross sections; as discussed 
in sect. 2, the generation of the hard process in MC@NLO has nothing to do with the 
analogous generation occurring in the MC. On the other hand, the factor of a s which 
explicitly appears in eqs. (5.6)-(5.18) is due to parton showering, and its argument 
should in principle be chosen according to eq. (A. 4). However, the choice of this 
argument is strictly speaking a matter beyond NLO. This is discussed more fully in 
I, sect. B.3. As done there, for simplicity we choose the same scale as that used for 
NLO terms. 

The analogous results for qq, qg, gq, and gq can easily be obtained from the 
equations above. The situation is summarized in table 2. 

These results allow one to obtain the MC subtraction terms c?S a f,| MC entering the 
MC@NLO definition, eq. (2.1), as explained in app. B. 

6. Implementation of MC@NLO 

The practical implementation of MC@NLO for heavy flavour production proceeds in 
a similar way to that for gauge boson pair production, described in sect. 4.5 of I. The 
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ab 


qq -> QQ 


qq -> QQ 


00 -> QQ 


qq 


±(t), Q(t), Q(t) 






qq 




±(w), Q(u), Q(m) 




qg 


-(*) 

V / 




+(s, £, u) 

1 \ 7 7 / 


qg 




-(«) 


+(s,*,u) 


gq 




+(«) 


-(s,t,u) 


gq 


+(*) 




-(s,f,u) 


gg 






±(s, Q(t,u), Q(t,u) 



Table 2: Short-distance contributions to MC subtraction terms, from Born processes 
qq — > QQ 5 99 — ^ QQi an d 55 — ^ QQ- Each entry lists the emitting legs (+, -, Q, Q); for 
each emitting leg, we report in parentheses the different contributions I, according to the 
possible colour flows (corresponding to Eq = \T\/2). 

integrals necessary to determine the required numbers of HI (2 — > 3 parton) and § 
(2^2 parton) configurations are computed according to eq. (1.4.23) and eq. (1.4.24), 
together with the corresponding equations in which the integrands are taken in ab- 
solute value. Configurations are generated according to the MC-subtracted weight 
distributions in eq. (2.1), and then unweighted with weight ±1 according to the sign 
of the weight, using the SPRING-BASES package [8]. Weighted events could also be 
generated, but this option is not implemented at the moment. 

The selected configurations and their weights are written on a file which is input 
to the HERWIG event generator. HERWIG then performs the parton showering and 
hadronization as explained in app. A. 2 and A. 3 respectively. The generated events 
with weight +1 are to be treated as real events for the purposes of histogramming 
and/or detector simulation, whereas those with weight -1 (a fraction of the order of 
10% for top and 20% for bottom) are "counter-events", which must be treated as 
real events in detector simulation but have to be subtracted from histogram bins to 
which they contribute. 

The only significant differences from the implementation of MC@NLO described 
in I are that a colour flow must be assigned to each parton configuration and that 
this flow, together with the event weight and the parton identities and momenta, are 
handled by the "Les Houches" generic user process interface [9] . Details of these two 
new features are given in the following subsections. 

6.1 Colour flow assignment 

The assignment of colour flows for § configurations is done according to the HERWIG 
prescription explained in app. A. 1.3: where the colour flow is ambiguous, it is assigned 
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Figure 4: Basic colour flow configuration involving three gluons and a heavy quark- 
antiquark pair. 






— > — 


1 


\ 

— < — 



Figure 5: Planar graphs contributing to the colour flow configuration of fig. 4. 



according to the N —>■ oo limit. The same prescription is extended to H configurations 
as follows. 

Consider first the process gg — > QQg. In the large-iV limit, the basic planar 
colour flow amplitude involving three gluons and a heavy quark-antiquark pair is 
depicted in fig. 4. Double-directed lines are gluons, and single-directed lines are 
quarks. The momenta 1%, I2 and Z3 are all defined to be outgoing. In the notation of 
sect. 3.3, the six independent colour-ordered amplitudes for the gg — > QQg process 
are obtained by replacing l lf l 2 , I3 — > —pi, —P2, k in all possible ways. 

The Feynman graphs that allow for the colour pattern of fig. 4 are shown in 
fig. 5. The amplitudes are easily computed from the large- N limit Feynman rules, 
which amount to the following graphical prescriptions for the vertices: 



qq gluon: 



1 

7! 



(6.1) 



three gluons: —= 
v2 



[g^{p? -p?) + g» m (rf -p^ + g^ip? -pT )] 
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Figure 6: Basic colour flow configuration involving a gluon, a heavy quark-antiquark pair 
and a light quark-antiquark pair. 



fourgluons: y (g^g^ 4 + g MW / 4 " 1 - 2g^ 3 g^ 4 ) , (6.2) 



where all momenta and indices have to be assigned counterclockwise, and momenta 
are all incoming. 
We define 



2 

(6.3) 



spin, col 

the square of the sum of the amplitudes in fig. 5 in the large- N limit, summed over 
spin and colours. The ordering of the momenta in the arguments of corresponds 
to the ordering of the colour connection, starting from the outgoing heavy quark. 

There are two basic colour orderings for processes involving a gluon, a light 
quark-antiquark pair, and a heavy quark-antiquark pair. They are depicted in fig. 6. 
Here we define two squared amplitudes 

f^(h, k 2 , l g , l q , l q ), f^fa, k 2 , l g , l q , l q ), (6.4) 

corresponding to colour flows (a) and (b) of fig. 6, respectively. The two squared 
amplitudes are related by charge conjugation 

f^(k 1: k 2 , l g , l q , l q ) = f^{k 2) k u l g , l q , l q ) . (6.5) 

Notice that in our notation we call q and q the light flavour lines corresponding to 
an outgoing quark or antiquark respectively. The corresponding colour flow squared 
amplitudes for qq — > QQg are obtained from the amplitudes of eq. (6.4) by the 
replacement l q = —pi, l q = —p 2 and l g = k. The squared amplitudes for the process 
qg — * QQq are instead obtained with the replacement l q = —pi, l g = —p 2 , l q = k. 

The correctness of our calculation of the colour flow amplitudes was checked by 
comparing the sum of all colour squared amplitudes with the known full squared 
amplitude, for very large values of N. 
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As an example, consider the assignment of colour flow for the process qg — > QQq. 
The configuration of momenta pi,P2, ki, k 2 , k is chosen according to the full (N = 3) 
expressions with modified subtraction as explained in sect. 2. Then the squared 
amplitudes for the large- iV colour flows in fig. 6 are computed: 

fa = f QmQ {K h, -P2, k, -Pi) , f b = f QmQ {k 2 , k h -P2, - Pl , k) . (6.6) 
The colour flows (a) and (b) are then assigned with probabilities 

J a + Jb Ja + Jb 

6.2 Interface to HERWIG MC 

The colour flow selected as described above is encoded in a single integer IC which 
is written on a file together with the event weight and the parton momenta and 
identities. The possible values of the colour flow code IC and their meanings are 
given in app. C. 

The file of parton configurations is read by HERWIG as input to the Les Houches 
generic user process interface [9]. A negative value of the process code IPROC signals 
to HERWIG that, instead of generating a partonic hard subprocess as outlined in 
app. A.l, it should load and use subprocess information in the Les Houches common 
blocks HEPRUP and HEPEUP. 

First a file header is read by the interface subroutine UPINIT, which checks that 
the file has been generated with parameter values consistent with those set in this 
HERWIG run and initialises the run common block HEPRUP. The parton configurations 
are then read sequentially by the interface subroutine UPEVNT, which loads the event 
common block HEPEUP with all information necessary to generate an event, including 
the interpretation of the colour flow code IC as outlined above. 



7. Results on top quark production 

In this section we present results for top quark production, obtained with the MC@NLO 
implementation 4 described above. We do not attempt here to give a phenomenolog- 
ical description that corresponds to a specific experimental configuration; rather, we 
wish to show the differences between MC@NLO, standard HERWIG MC, and NLO 
results. For this reason, we present distributions obtained either by integrating over 
the whole phase space, or else with cuts on the heavy quark variables, rather than on 
those of their decay products. We first show results for the LHC, i.e. for pp collisions 
at v^S* = 14 TeV. We set the top mass to m — 173 GeV, and the renormalization 
and factorization scales equal to the top transverse mass, ^Jm? + p\. We adopt the 

4 A public version of the program is available on the MC@NLO web page 
http : //www . hep . phy . cam . ac . uk/theory/webber/MCatNLO/. 
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low-a s set of MRST99 parton densities [10], since this set has a A QCD value which is 
rather close to that used as HERWIG default. This value of A QCD (Alf = 164 MeV) is 
used for all our MC@NLO, MC, and NLO runs. In the case of standard HERWIG MC 
runs, we give each event (we generate unweighted events) a weight equal to crtot/N to t, 
with N to t the total number of events generated. For fully exclusive distributions with 
no cuts applied, this is equivalent to normalizing HERWIG results to the total NLO 
rate a tot (with our choice of parameters, we obtain <j tot = 6.668 pb at the Tevatron, 
and a tot = 736.6 pb at the LHC). All the MC@NLO and MC results (but not, of 
course, the NLO ones) include the hadronization of the partons in the final state. 
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Figure 7: MC@NLO (solid), HERWIG (dashed) and NLO (dotted) results for the rapidity 
(left panel) and the transverse momentum (right panel) of the top quark at the LHC, with 
acceptance cuts. HERWIG results have been normalized as explained in the text. 



We present in fig. 7 the rapidity (left panel) and the transverse momentum (right 
panel) distributions of the top quark. The rapidity (transverse momentum) result has 
been obtained after applying the cut > 20 GeV (\y (t) \ < 1). The solid, dashed and 
dotted histograms show the MC@NLO, MC and NLO results, respectively. These 
distributions are fairly inclusive, and we expect them to be reliably predicted by NLO 
QCD in a wide range. From the figure, we see that the NLO and MC@NLO results 
are extremely close to each other in the whole ranges considered. HERWIG results 
are also fairly close to the NLO and MC@NLO ones, giving only a slightly broader 
rapidity distribution. The same pattern can be found for the invariant mass of the ti 
pair. We conclude that, for these kinds of observables, NNLO effects are very small, 
and NLO, MC and MC@NLO are almost equivalent. This also implies that any 
possible reshuffling of the momenta, due to the hadronization phase in MC@NLO, 
has negligible impact on the t and i. The right panel of fig. 7 also clearly shows 
a characteristic feature of the comparisons between MC@NLO and NLO results, 
namely that the former are numerically more stable than the latter. This feature, 
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which is even more evident in the case of b production, is due to the fact that, as 
eq. (2.1) documents, in MC@NLO all cancellations between large numbers occur at 
the level of short-distance cross sections, rather than in histograms as in the case of 
NLO computations. 

We now turn to the case of more exclusive quantities, such as correlations between 
t and t variables. In fig. 8 we present the modulus of the vector sum of the transverse 
momenta of the t and t, which we denote by p^ . NLO computations cannot predict 
this observable in the region ~ 0, because of a logarithmic divergence for ~~ > 
0; on the other hand, NLO is expected to give reliable predictions at large . The 
MC behaves in the opposite way; thanks to the cascade emission of soft and collinear 
partons, it can effectively resum the distribution around p? — 0. However, its 
results are not reliable in the large-p^ region, which is mainly populated by events 
in which a very hard parton recoils against the tt pair. 




log 10 (p ( f/GeV) 3 1 log 10 (p ( ^/GeV) 



Figure 8: As in fig. 7, for the transverse momentum of the tt pair, without (left panel) 
and with (right panel) acceptance cuts. 

The complementary behaviour of the NLO and MC approaches can be seen 
clearly in fig. 8, regardless of the cuts on the rapidities and transverse momenta 
of the heavy quarks. In the tail of the p^ t} distribution, the NLO cross section is 
much larger than the MC one, simply because hard emissions are correctly treated 
only in the former. The presence of the dead zones shown in figs. 2 and 3 makes it 
very difficult to generate a very hard parton recoiling against the tt pair in the MC, 
which therefore gives rates much below the NLO result in this region. For p£ ~~ * 0, 
the difference between the two histograms shows the effect of all-order resummation; 
clearly, no meaningful comparison between NLO and data can be attempted in this 
region. It is therefore reassuring that the MC@NLO result interpolates the MC and 
NLO results smoothly. In the small-p?*' region, the shape of the MC@NLO curve 
is identical to that of the MC result. This is evidence of the fact that MC and 
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MC@NLO resum large logarithms at the same level of accuracy, as argued in I. 
When p$ grows large, the MC@NLO tends to the NLO result, as expected. Again, 
hadronization has no significant impact on the ti system. 




Figure 9: As in fig. 8, for the difference in the azimuthal angles of the t and t. 



In fig. 9 we present the distribution of the difference between the azimuthal 
scattering angles (i.e., those in the plane transverse to the beam direction) of the t 
and i, which we denote by A0 (tt) . This distribution cannot be reliably predicted by 
fixed-order QCD computations in the region A0 (tt) ~ it; in fact, the NLO prediction 
diverges logarithmically for A0 (tt) — > n. This can be seen in the insets of the plots, 
where the cross section has been plotted versus (ix — A0 (tt) )/ tc on a logarithmic scale, 
in order to visually enhance the region A0 (tt) ~ ti. We can see that in this region the 
MC@NLO and MC results have identical shapes, as in the case of the observable p^ 
near zero discussed above. The other end of the spectrum, i.e. the tail A0 (tt) ~ 0, 
is populated by configurations in which a hard jet recoils against the ti pair, but 
also by configurations in which the t and i have small transverse momenta. The 
NLO calculation gives a good description of the former region, but not the latter, 
while the MC can treat reliably the latter region, but not the former. MC@NLO, on 
the other hand, is expected to handle both regions reliably, while still avoiding any 
double counting. We also notice that, when cuts on the transverse momenta of the 
heavy quarks are applied, the contribution from multiple soft or collinear emissions 
becomes less important, and so the MC does less well at A0 (tt) ~ 0, as may be seen 
in the right-hand panel of fig. 9. 

Corresponding results for Tevatron Run II {pp at \/S = 2 TeV) are shown in 
figs. 10-12. An interesting new feature of the top quark rapidity distribution is its 
forward-backward asymmetry, which cannot appear in pp collisions and is also absent 
in pp at the Born level, and hence also in the MC result. To document this, we plot 
in the left panel of fig. 10 the rapidity asymmetry rather than the rapidity itself; cuts 
Pt ,Pt > 20 GeV have been imposed. We see a fair agreement between MC@NLO 



28 



and NLO, whereas HERWIG predicts an asymmetry compatible with zero. In the 




log 10 (pV/GeV) 



Figure 10: MC@NLO (solid), HERWIG (dashed) and NLO (dotted) results for the ra- 
pidity asymmetry (left panel) and the transverse momentum (right panel) of the top quark 
at the Tevatron. HERWIG results have been normalized as explained in the text. 



right panel of the same figure, the top quark transverse momentum (with a \y (t) \ < 1 
cut) is shown. As in the case of LHC, we see no substantial differences between NLO, 
MC@NLO and HERWIG. 
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Figure 11: As in fig. 10, for the transverse momentum of the tt pair, without (left panel) 
and with (right panel) acceptance cuts. 



The picture for (fig. 11) is broadly similar to that at the LHC, except for 
the reduced tail at high p^ l) : the MC@NLO prediction makes a smooth transition 
from the NLO to the resummed MC form as p 1 ^ decreases. 

The situation for A0 (tt) , on the other hand, is slightly different. At Tevatron 
energies the influence of hard emissions is not so strong as at the LHC. Consequently 
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(see fig. 12) the MC and MC@NLO predictions coincide quite closely over the whole 
range of A0 (tt) , whereas the NLO remains lower, regardless of whether or not cuts 
are applied. 



8. Results on bottom quark production 

In the case of bottom production, many complications arise that are not present 
in top production. In a standard Monte Carlo, all of the three mechanisms dis- 
cussed in sect. 3.1, namely flavour creation (FCR), flavour excitation (FEX), and 
gluon splitting (GSP), need to be considered; they have rather different kinematic 
signatures, and they are dominant in different regions of the phase space. From the 
point of view of perturbative computations, bottom cross sections are characterized 
by a fairly large value of the coupling constant, which implies sizable K factors; also, 
the importance of large logarithmic terms arising at all orders is manifest in many 
observables, and suitable resummations are often necessary for sensible comparisons 
with data. 

As in the case of top production, we shall not attempt here to discuss the phe- 
nomenological implications of our findings. We shall rather emphasize the kind of 
problems encountered in 6-physics simulations with standard MC's, and the way in 
which MC@NLO solves some of them. Comparisons between MC@NLO and NLO 
results will also be presented here. We shall not discuss the well-known limitations of 
fixed-order computations, but refer the reader to ref. [11] for further details. Unlike 
the case of top production, in b physics we can compare MC@NLO results with MC 
and NLO ones either for b quarks or for 6-flavoured hadrons. The former option 
is clearly preferable, if we aim at understanding the extent of the validity of the 
MC@NLO approach, and the improvements with respect to traditional formalisms. 
Technically, in MC@NLO and in HERWIG we define the 6-quark level as the stage 
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which comes immediately before the gluon-splitting phase, 5 corresponding to HER- 
WIG status code ISTHEP=2. 

The results presented in this section have been obtained at the Tevatron Run II 
(pp at v^S* = 2 TeV), for a bottom mass of 5 GeV. The other parameters have been 
chosen as in sect. 7. In order to simplify the analysis procedure, pair observables 
are defined by considering all possible pairs in the event, regardless of their charge. 
Thus, in HERWIG and MC@NLO bb, bb, and bb pairs are treated on an equal footing 
(at the NLO, one has just one bb pair), and will collectively be denoted as bb pairs. 
From the practical point of view this detail is almost irrelevant, since we find that 
the probability of having more than two fe's in an event with at least two fe's is of the 
order of 0.1%. 

8.1 6-production issues in HERWIG 

We start by discussing the problems arising in the simulation of b production with 
HERWIG. We stress that similar problems are present in any standard parton shower 
MC. As we discussed in sect. 3.1, FEX and GSP contributions are considered in 
heavy flavour production simply because FCR alone is not capable of describing the 
kinematics of observed events. It should be noted that FEX and GSP are somewhat 
anomalous from the point of view of MC's, since usually the simulation relevant to a 
given hard system involves the production of such a system at the level of hard process 
generation. This fact has profound consequences: MC's cannot simulate small-p T 
production of heavy quarks, since FEX and GSP matrix elements are diverging for 
p T — > 0. This poses a practical problem, which is easily circumvented by cutting 
off the matrix element divergencies. In HERWIG, this is achieved by requiring the 
transverse momenta of the primary partons to be larger than a given quantity, called 
PTMIN. In addition to this, HERWIG has an effective cutoff at the level of hard matrix 
elements for FEX processes, which prevents the generation of primary partons with 
p T = even if PTMIN = (see app. A. 1.2). In b production, GSP processes also have 
an effective cutoff, but of a different nature. The probability of getting a showering 
scale large enough to produce a bb pair vanishes as p T — > in the hard process. Still, 
this doesn't allow one to set PTMIN = in GSP, since the hard process is generated 
independently of the shower. 

Although t-channel singularities are cut off by PTMIN, a problem of principle 
remains: if one interprets the output of a given showering event as a Feynman 
diagram, one obtains a contribution which in perturbation theory is only relevant in 
the large-momentum regime p T 3> m. Thus, strictly speaking one should run FEX 
and GSP, and keep only those events with a large-p T heavy quark. It is customary to 
ignore this problem, and to keep all the events generated. The results are in general 

5 This gluon splitting is non-perturbative, and it has nothing to do with the perturbative branch- 
ing of a gluon which takes place during a shower; see app. A. 3.1 for more details. 
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biased by PTMIN, and it is therefore necessary to insure that this bias does not affect 
the predictions in the kinematical regions of interest. 

We have studied the bias due to the choice of PTMIN by considering the GSP 
contribution. An analogous study can be done for the FEX contribution, but in 
this case the presence of an effective cutoff at the level of hard matrix elements 
complicates the discussion unnecessarily. We have considered the inclusive b cross 
section, requiring 

(8.1) 



p£> > 5 GeV , 



y (h) \ < 1 



and the pair cross section, requiring 



p^\p? > 5 GeV : 
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The results are presented in fig. 13; open points are HERWIG predictions for the 
corresponding PTMIN value, and the solid lines are there just to guide the eye. The 
dotted line is the weighted average of the results for the pair cross sections obtained 
at PTMIN = 3, 5, and 6 GeV. The lower panel of fig. 13 has an easy interpretation: for 
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Figure 13: b (upper panel) and bb (lower panel) GSP cross sections within the cuts given 
in eqs. (8.1) and (8.2) respectively, as predicted by HERWIG for various choices of PTMIN. 



PTMIN values smaller than 6 GeV, HERWIG predictions are independent of PTMIN, 
i.e., they are not biased. Clearly, these values of PTMIN are correlated to the choices 
of cuts made in eq. (8.2) and, to a smaller extent, to the fact that the observable 
chosen is fully inclusive. It seems however safe to choose PTMIN = 5 GeV for studying 
any type of pair correlations with the cuts of eq. (8.2). 
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The upper panel of fig. 13 displays a less pleasant behaviour; no range in PTMIN 
can be found where the inclusive-6 cross section is independent of PTMIN. This 
happens because, at small but finite p T of the primary partons, HERWIG can choose 
a colour flow according to which the evolution scale is almost equal to s. Thus, a 
gluon acquires a large enough virtuality to split into a bb pair. When boosted to 
the lab frame, one of the 6's can have a transverse momentum exceeding the p T cut 
(the probability for this to happen is small, and thus the probability of getting both 
6's above the cut is negligible, which explains the difference between the two panels 
of fig. 13). As already discussed before, the selection of such a colour flow is less 
and less probable with decreasing p T . Thus, there must exist a PTMIN which returns 
unbiased single-inclusive cross sections. We didn't try to find such a PTMIN value 
here, since for single-inclusive distributions it is more sensible to compare MC@NLO 
results to NLO ones. However, our exercise proves that it is very time-consuming to 
get unbiased predictions with HERWIG: the efficiency for generating events passing 
the cuts of eqs. (8.1) and (8.2) is very rapidly decreasing with decreasing PTMIN, 
because of the divergence of the matrix elements at PTMIN = 0. 

To document the relative importance of the mechanisms contributing to HER- 
WIG predictions, we show in fig. 14 the results for the azimuthal distance A^ (bb) 
between the b and b, and the transverse momentum p^ b) of the bb pair. The dashed, 
dotted, and dot-dashed histograms are the FCR, FEX, and GSP contributions re- 
spectively, whereas the solid histogram is the sum of the three. It is apparent that 
FCR is important only for those kinematic configurations which are almost 2 —>■ 2, 
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Figure 14: Azimuthal bb distance (left panel) and transverse momentum of the bb pair 
(right panel) as predicted by HERWIG (solid); the contributions of FCR (dashed), FEX 
(dotted), and GSP (dot-dashed) are also separately presented. 

namely A0 (bb) ~ n and p^ h) ~ 0; elsewhere, FEX and GSP contributions cannot 
be neglected. This implies that 6-physics simulations in standard MC's are always 
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computing intensive. We performed our GSP runs setting PTMIN = 5 GeV; we found 
that the fraction of events with at least a bb pair is about 3.7- 10~ 3 ; when the cuts of 
eq. (8.2) are applied, this fraction is reduced to 1.6-10 -4 . The dot-dashed histograms 
in fig. 14 contain only 1900 events, obtained by generating 1.2-10 7 GSP events with 
HERWIG. Although techniques are known to increase the efficiency, and the figures 
quoted above depend on the particular MC used, it is at present unknown how to 
generate GSP bb events with high efficiency 

8.2 6-production issues in MC@NLO 

In standard shower Monte Carlo programs, the bottom flavour is included in the 
evolution of the initial-state and final-state showers, consistently with the fact that 
FEX and GSP processes need also to be considered. On the other hand, as discussed 
in sect. 3.1 the implementation of MC@NLO is based on FCR processes only. Thus, 
if we implemented bottom production in exactly the same way as we did with top, 
we would have to switch off the bottom flavour in all places where it appears in 
the shower Monte Carlo. This would be particularly problematic for the initial-state 
shower. Since we normally use five-flavour parton densities, switching off the bottom 
flavour in backward evolution may lead to inconsistencies. 

This problem is not peculiar to MC@NLO implementation. The original NLO 
calculation of heavy flavour production was carried out in the decoupling scheme of 
ref. [12], and thus should be used in conjunction with a 4-flavour coupling constant 
and 4-flavour parton densities. In practice, five-flavour parton densities were always 
used, since, as pointed out in eq. (3.11) of ref. [13], this is correct (up to a numerically 
small effect in the gluon parton density) as long as one neglects parton densities with 
incoming heavy quarks. 

In ref. [6], the exact prescription for a change of scheme in the heavy flavour 
production formulae (in order to go from the decoupling scheme of ref. [12] to the 
full MS scheme with 5 flavours) was given. It is summarized as follows: 

• Use five-flavour parton densities and strong coupling constant. 

• Ignore the b and b flavours in the parton densities. 

• Add a term —a s ^-\og ^ (where n R is the renormalization scale) to the 
qq channel cross section. 

• Add a term — a s ^Mog^ afg (where /i P is the factorization scale) to the gg 
channel cross section. 

This change of scheme was implemented in the MC@NLO code. The corrections in 
the last two items above, although necessary to maintain formal correctness at the 
NLO level, have very small numerical impact. The implementation of the MC@NLO 
for bottom is then identical to that for top, except for the corrections listed above. 
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Unlike the case of top production, typical studies of b production at hadron 
colliders are carried out in a relatively large transverse momentum regime, where 
the bottom flavour behaves in part as a light flavour, and the resummation to all 
orders of L = log p T /m terms may therefore be necessary in order to obtain sensi- 
ble predictions. In perturbation theory, and for single-inclusive observables, such a 
resummation is carried out to NLL accuracy by considering all 2 — > 2 and 2 — > 3 
partonic processes (i.e., including those with one or no b in the final state), and con- 
volving them with NLL perturbative fragmentation functions [14], which describe 
the evolution of a light parton or a b quark into a b quark. The resummed formulae 
can be consistently combined [6] with the exact 0(a^) formulae, in order to obtain 
reliable predictions for any p T value. 

This procedure has not been extended yet to more exclusive observables. MC@NLO 
would seem a natural way to achieve this goal, at least at the leading logarithmic ac- 
curacy; this would imply implementing, besides the exact 0{a\) FCR formulae, also 
the O(otf) ones for light-parton scatterings. The role of the perturbative fragmen- 
tation functions would then be played by the showers. The presence of light-parton 
processes would be taken into proper account in the definition of the MC subtrac- 
tion terms, in order to avoid double counting. Although this task is not totally out 
of reach, it is very difficult to achieve it in practice, since it would lead to a great 
increase in complexity in comparison to the top case. 

In the present work we shall restrict ourselves to implementing bottom produc- 
tion with the same accuracy as standard NLO calculations, using the five-flavour 
strong coupling constant and parton densities in the scheme described above. This 
allows us to include some leading-log effects (i.e., terms like al(a s L) k ), but not all 
of them. More precisely, single-log terms (of order a^L) are included exactly in our 
calculation. The terms of order a^(a s L) k that are included for any k are those that 
result from multiple final-state radiation from the heavy quark lines, one example of 
which is illustrated in fig. 15. Observe that, although the hardest emission is exactly 



Figure 15: Enhanced terms present in the MC@NLO implementation, accurate to the 
leading-logarithmic level in the MC@NLO approach. The shaded area represents an almost 
collinear branching due to parton showers. 

described at the NLO level, subsequent gluon emissions have only a leading-log ac- 
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curacy, since going to next-to-leading log accuracy would require the implementation 
of emission kernels (i.e., of showers) at comparable accuracy. 

There are also other effects of order ag(a s L) k that are automatically included in 
our MC@NLO implementation, which are due to branchings that involve b quarks. 
One such contribution arises when, by backward evolution, a gluon fusion FCR 
process is connected to a b quark line via the initial-state shower. A second type 
of contribution arises when a final-state gluon, emerging from a FCR hard process, 
splits into a bb pair. This splitting process, together with the associated virtual 
process, gives rise to the correction to the running of a s due to the bottom flavour. 
The consistency between the shower evolution including b quarks, and the NLO FCR 
cross sections, is in fact what dictates the use of the five-flavour scheme for a s and 
for the parton densities. 

However, it is clearly impossible to include all effects of order a>l(a s L) k starting 
from FCR processes, even with a five-flavour scheme. Most noticeably, the al(a s L) 2 
GSP and FEX contributions depicted in fig. 1 are not included in MC@NLO. 

In the case of final states containing two or more heavy flavour pairs, the level 
of precision of the MC@NLO treatment is unclear because one of the pairs (the one 
produced at the primary NLO vertex) is described differently from those produced 
by parton showering, as illustrated in fig. 16. As a consequence, the weight factors 
in the prediction of the bottom pair multiplicity, for example, may not be correct. 
We postpone discussion of this point to a later date, since from a practical viewpoint 
multiple pair production is a rare phenomenon (0.1% of single-pair production) which 
has a negligible effect on the plots shown here. 




Figure 16: Heavy flavour production by gluon splitting, followed by a further splitting 
process (by showering) of the recoiling gluon, represented by a shaded area. 



8.3 Pair correlations 

In this section, we present predictions for bb pair observables. We start by considering 
the transverse momentum of the pair, p^ b) . The MC@NLO (solid) and NLO (dotted) 
results are shown in fig. 17; in the left panel no cuts have been applied, whereas the 
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right panel includes the effect of the cuts given in eq. (8.2). Regardless of the presence 
of the cuts, the two plots display the same pattern as the analogous plots for top 
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Figure 17: MC@NLO (solid), HERWIG (dashed), and NLO (dotted) results for the pair 
transverse momentum, without (left panel) and with (right panel) the cuts of eq. (8.2). 



production (figs. 8 and 11): the NLO distributions diverge for p ( ^ h) — > 0, whereas 
MC@NLO results have a regular behaviour. On the other hand, in the large-px b) 
region, dominated by hard-parton emission, MC@NLO and NLO coincide in shape 
and normalization. 

When the cuts of eq. (8.2) are applied, we can also consider the HERWIG result 
(dashed histogram), which has been already shown in fig. 14. Unlike the case of 
top production, here we take HERWIG result with its normalization; a rescaling by 
the K factor would be more appropriate only if FCR process alone were included. 
The agreement between MC@NLO and HERWIG is remarkable. In the low-px b) 
region, we can see that the shape of the MC@NLO result is basically the same as 
that of HERWIG, analogously to what we have seen in the case of top production. 
The comparison at large p ( ^ b) is hampered by the poor statistics of the HERWIG 
result, since this region receives its main contribution from the GSP process (see 
fig. 14); however, it appears that MC@NLO and HERWIG agree well. This is not 
surprising, since a bb pair produced through the GSP mechanism mainly recoils 
against a hard gluon, exactly as in the NLO matrix elements which are implemented 
in MC@NLO. The inability of the MC to produce hard emissions is only evident in 
the FCR contribution (dashed histogram in fig. 14), which in fact has a much softer 
behaviour than GSP or MC@NLO. 

We now turn to the case of the azimuthal distance, A0 (bb) . As discussed in 
sect. 7, this observable receives non-negligible contributions to the tail A0 (bb) ~ 
from both hard emissions and multiple soft or collinear emissions. The presence 
of large logarithms in the perturbative expansion is manifest at the NLO, since the 
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prediction is strongly peaked towards A0 (bb) — ► 7r; the correct result for the total rate 
implies that the cross section at A0 (bb) = tt is negative. This behaviour is apparent in 
the dotted histograms in fig. 18. The shapes of the MC@NLO results are similar to 
the NLO ones in the small-A0 (bb) region, but do strongly differ elsewhere. One could 
perhaps expect that, upon applying the cuts of eq. (8.2), MC@NLO and NLO would 
be closer in normalization for A0 (bb) — > 0; this doesn't happen because MC@NLO and 
NLO have non-negligible differences in the intermediate-p T region, which dominates 
the cross section when the cuts of eq. (8.2) are applied. This fact is not accidental, 
and will be discussed in sect. 8.4. 

As in the case of the transverse momentum of the pair, when the cuts are ap- 
plied we can also consider the HERWIG prediction (dashed histogram in the right 
panel of fig. 18). In this case, too, the agreement between MC@NLO and HERWIG 
is remarkable, and emphasises again the importance of the GSP contribution (see 
fig. 14). MC@NLO appears to be somewhat more peaked than HERWIG towards 

A0 (bb) -> TT. 



8.3.1 Impact of initial-state radiation 

The results presented above imply that multiple radiation is a crucial effect in b 
production. Lacking an N fc LO computation, a typical way of estimating the effects of 
multiple initial-state emissions on heavy-quark distributions is to supplement an NLO 
calculation with an intrinsic transverse momentum for the incoming partons (denoted 
as NLO+A; T -kick hereafter). We stress that this procedure is ill-defined, since there is 
no way of avoiding double counting. There are however cases in which it is justified 
from the phenomenological point of view, allowing a much better description of the 
data than NLO predictions alone (see ref. [11] for a discussion of the implementation 
of /c T -kick in heavy flavour production, and its implications). 
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Figure 19: A0 (b6) (left panel) and p? b) (right panel) predicted by MC@NLO (solid), 
NLO+fcr-kick (dashed), and NLO (dotted). The cuts of eq. (8.2) have been applied. 



We define our /c T -kick by assuming a gaussian transverse momentum distribution 
for the incoming partons; we fix the free parameter of this distribution by requiring 
that the average p { ^ b) predicted by NLO+/c T -kick be equal to that obtained from 
MC@NLO. In this way, the average transverse momentum of the incoming partons 
turns out to be (k T ) ~ 4 GeV. With this choice, the NLO+/c T -kick results for the 
£>T b) and A0 (bb) distributions are seen to agree reasonably with the corresponding 
MC@NLO results (see fig. 19). This confirms the importance of multiple parton 
radiation, on top of the possible hard emission present at the NLO level, and shows 
that the vast majority of the effect is due to emissions from initial-state partons. 

The reader is urged not to take NLO+/c T -kick results too seriously, since they 
are based on a model with little theoretical justification. In fact, an average in- 
trinsic transverse momentum of 4 GeV is too large to be considered a typical non- 
perturbative effect. On the contrary, in the MC@NLO implementation this effect 
has a purely perturbative origin, being due to multiple initial-state emissions. 

8.4 Single-inclusive observables 

We finally turn to the case of single-inclusive 6-quark distributions. In this section, we 
shall not consider HERWIG results, because of the findings of sect. 8.1; on the other 
hand, we do consider the effect of the NLL resummation of large logs L = logp T /m, 
as discussed in sect. 8.2. 

We start with the rapidity of the b quark, presented in fig. 20 without and with 
a transverse momentum cut > 5 GeV. Regardless of the presence of this cut, 
MC@NLO (solid) and NLO (dotted) results agree well. This is to be expected, since 
for such an observable NLO results are in general reliable. 
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Figure 20: MC@NLO (solid) and NLO (dotted) results for single-inclusive b rapidity 
spectrum, without (left panel) and with (right panel) a p T cut. 
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Figure 21: MC@NLO (solid) and NLO (dotted) results for single-inclusive b p T spectrum, 
without (left panel) and with (right panel) a rapidity cut. In the latter case, the FONLL [6] 
result (solid line) is also shown. 



The differences between MC@NLO and NLO results are larger in the case in 
which a p T cut is applied. The reason becomes clear if one plots the p T spectrum, 
shown in fig. 21. Regardless of the presence of the rapidity cut, the MC@NLO 
predictions (solid) are above the NLO ones (dotted) in the intermediate-p T range. 
When we apply a rapidity cut, we can also consider FONLL predictions [6], based on 
a formalism that includes not only NLO, 0(al) terms, but also log-enhanced terms 
of order a\{a s \ogp T / m) k (LL) and a\{a s \ogp T /m) k (NLL). The FONLL result is 
shown as a solid line in the right panel of fig. 21. 

Similarly to MC@NLO, the FONLL result is slightly above the NLO calculation 
in the p T region from about 5 to 50 GeV. In the FONLL calculation, this small 
excess in the intermediate region was attributed to the inclusion in the calculation of 
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higher-order corrections to light parton production (in particular gg — > gg) followed 
by gluon fragmentation into a heavy quark (see ref. [6] for details). Neither the 
NLO, nor the MC@NLO calculations include these corrections. On the other hand, 
in the MC@NLO implementation the enhancement is very likely due to the transverse 
momentum boost given to the hard process by the initial state showers. In order to 
verify this possibility, we use again the NLO+/c T -kick approach used in sect. 8.3.1, 
with the same parameter setting adopted there. The result is presented in fig. 22. 
The MC@NLO and NLO+/c T -kick results are again in fair agreement. As in the 
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Figure 22: Comparison of the 6-quark p T spectrum obtained with MC@NLO (solid), 
NLO+£; T -kick (dashed), and NLO (dotted). 

previous section, we interpret this fact as a hint of the relevance of initial-state 
multiple emissions in b production at the Tevatron. Such effects are not included in 
FONLL. We thus tentatively speculate that a formalism able to take into account 
all these effects (i.e. FONLL effects plus multiple soft radiation from the initial- 
state partons) may result in further enhancement of the intermediate-p T region with 
respect to both MC@NLO and FONLL results. 

9. Conclusions and future prospects 

In this paper we have applied the MC@NLO method to the process of heavy quark 
production at hadron colliders. The method was developed in ref. [1] and applied 
there to the process of gauge boson pair production. For heavy quark production, 
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the basic formalism remains the same; it was only necessary to compute the relevant 
MC subtraction terms and associated colour flows. However, this proved a good 
deal more complicated than for gauge boson pair production, owing to the presence 
of coloured objects in the final state even at the Born level, giving rise to different 
patterns of momentum reshuffling in the MC and more alternative colour flows. In 
fact, the implementation for heavy quark production demonstrates that MC@NLO 
can be applied to arbitrarily complicated processes, in terms of kinematics, colour 
flows and large K factors. 

Results were presented on single-particle distributions and correlations at the 
heavy quark level, with parton showering and hadronization but showing the distri- 
butions of the heavy quarks before decay (in the case of top) or hadronization (in 
the case of bottom), to facilitate comparisons with purely perturbative predictions. 

In the case of top quark production, results appear to be fully under perturbative 
control at Tevatron and LHC energies. The MC@NLO results on single-particle 
distributions are generally close to those at NLO. A forward-backward asymmetry is 
seen in the top-quark rapidity distribution at the Tevatron, as expected in NLO but 
not in LO or in standard MCs which lack exact NLO corrections. In tt correlations, 
MC@NLO combines the more correct treatment of hard emissions in NLO with the 
MC resummation of soft and collinear contributions, giving a smooth distribution of 
Pt^ and an enhancement in the A0 (tt) distribution when the t and t are not back-to- 
back. The effects of hard emissions are naturally more visible at the LHC. 

Results on b production (presented for the Tevatron only) show more markedly 
the advantages of the MC@NLO approach compared with conventional MCs. The 
fact that NLO contributions are included eliminates the need for separate flavour 
creation (FCR), flavour excitation (FEX) and gluon splitting (GSP) contributions, 
which are required in standard MCs but are plagued with inefficiencies, ambiguities 
and cutoff dependences. In fact, we did not attempt to generate single-6 distributions 
with the standard HERWIG MC for comparison with MC@NLO and NLO, due to 
its low efficiency for the GSP contribution. The MC@NLO results, on the other 
hand, tend to be even more stable than pure NLO, since cancellations between large 
numbers occur at the matrix-element level, rather than in histograms. They do not 
depend on cutoffs or on b PDFs close to threshold (a source of uncontrolled errors), 
and thus one can sensibly predict cross sections all the way down to p T = 0. 

MC@NLO predictions for single 6-quark distributions at the Tevatron are quite 
similar to those at NLO, with some enhancement at intermediate transverse mo- 
menta, which can be interpreted as resulting from Sudakov effects due to multiple 
initial-state gluon emission. 

In the case of bb correlations at the Tevatron, we could obtain stable HERWIG 
MC predictions for reasonable values of the cutoff parameter PTMIN, provided the 
cuts of eq. (8.2) were imposed. With the resulting mixture of FCR, FEX and GSP 
contributions, the HERWIG results were then in quite good overall agreement with 
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those of MC@NLO. However, as stressed above, the arbitrariness and inefficiency of 
the standard MC prescription make the MC@NLO approach vastly preferable from 
all points of view. 

One feature of MC@NLO that might seem unattractive compared with standard 
MCs is the presence of a fraction (10-20%) of negative weights. As may be seen 
from the histograms presented here (all generated with unweighted events) this does 
not cause problems in practice. Most of the negative weights arise from the small-p T 
region where the NLO real emission contribution is not resolvable. It may be possible 
to reduce their contribution further by tuning the subtraction procedure. In spite of 
the presence of negative weights, MC@NLO is more efficient that standard MCs for 
b production, since it generates the GSP contribution with efficiency 1. 

We have not attempted any detailed phenomenological study or comparisons 
with experimental data in this paper, since our primary objectives were to present 
the MC@NLO method and compare it with others in the context of heavy quark 
production. An obvious next step would be to compare with Tevatron data on final- 
state properties in top and bottom production. In the case of top production, this 
will require the inclusion of decay correlations due to polarization of the t and t, 
which is straightforward in principle but not yet included in the NLO calculation 
that we have used. In the case of bottom production, careful attention must be paid 
to the hadronization model used to connect 5-quark and B-hadron distributions (see 
app. A. 3. 2 and ref. [15]). 

Other obvious future objectives are the application of MC@NLO to Higgs bo- 
son, single vector boson, and jet production. The only extra complication, compared 
with heavy quark production, arises in the latter case from collinear singularities 
due to emission from massless final-state partons. The necessary MC subtraction 
terms can be computed using the kinematics already presented in the present pa- 
per. The resulting terms will cancel the additional singularities and should provide 
more reliable predictions of jet production and fragmentation when combined with 
a subtraction-method NLO calculation of this process. 
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Appendices 



A. Heavy flavour production in HERWIG 

We describe here the main features of the HERWIG event generator that are relevant 
to heavy flavour production. We should distinguish between the standard HERWIG 
Monte Carlo program and the HERWIG implementation of parton showering and 
hadronization that forms part of MC@NLO. In the former case the program includes 
the generation of a hard subprocess as outlined in app. A.l. In the case of MC@NLO, 
only app. A. 2 and A. 3 are relevant, since the hard process configurations are read 
from the input file as explained in sect. 6.2. 

A.l Hard subprocess 

The primary heavy quark processes included in the standard HERWIG program are 
the 2 — ► 2 processes of flavour creation (FCR, eq. (3.1)), and flavour excitation (FEX, 
eq. (3.3) and charge-conjugate processes). However, the hard process configuration 
is always chosen as if it were FCR, i.e. the kinematics are always chosen for pip 2 — > 
ki, k 2 with pi 2 massless and /c 12 having the heavy quark mass m. 

A. 1.1 Kinematics 

The algorithm for selecting the 2^2 subprocess kinematics is as follows: 

1. Choose the heavy quark transverse mass m T = \/k 2 + m 2 from a power distri- 
bution da/dm T oc m~ p where p =PTP0W[4] 6 with PTMIN[10]< k T <PTMAX[ v / 5/2]. 

2. Choose the outgoing parton rapidities uniformly in YJMIN[-8]< 2/3 4 <YJMAX[8] 
(or in the kinematically allowed region if smaller). 

3. Compute the incoming parton light-cone momentum fractions according to 



x x = (e y3 + e V4 )m T /VS , x 2 = (e~ y3 + e^jm^VS . 



(A.l) 



The program also computes at this stage the cm. scattering angle 



cos# = (e ys - e y4 )m T / v / ey^(x 1 x 2 S - Am 2 ) . 



(A.2) 



4 



Compute the kinematic invariants according to 



s = X\X 2 S , 



t = -{l + e 



,V4-y3 



)m, 



2 



u = —s — t . 



(A.3) 



6 



Default values (in GeV units) are given in square brackets. 
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A. 1.2 Dynamics 

Having chosen a kinematic configuration, HERWIG calls for PDF evaluations at 
momentum fractions x 1 ,x 2 and hard process scale EMSCA given by 



V2stu , . . 

which is also used as the argument of a s in the calculation of hard process matrix 
elements. 

The program then cycles through all possible FCR and FEX subprocesses, com- 
puting the product of relevant PDFs times the on-shell differential cross section 
evaluated at (s, i, u) , times a Jacobean factor. The event weight EVWGT is the sum 
of all these contributions. In unweighted event generation, if this is greater than 
the maximum weight times a random number, the configuration is accepted. The 
program then multiplies EVWGT by another random number, and cycles through the 
FCR and FEX subprocesses again until the partial sum of contributions exceeds this 
value. The last subprocess selected is then loaded into the event record as an on-shell 
2^2 process, with the parton momenta reconstructed from the Xi,x 2 ,cos9 values 
computed previously. 

Note that: 

1. Values of (s, t, u) are generated only inside the FCR phase space, and so values 
accessible only in FEX are not explored. In particular the singular FEX point 
t — is never generated, even if one sets PTMIN=0. 

2. If the process selected is FEX instead of FCR, the reconstructed p T = p™ x of 
the heavy quark is not equal to the p T = p^ cn originally selected assuming FCR 
kinematics. An elementary calculation shows that 

y (s + m z ) [s — Am z ) 

Thus one always has p™ x > p^ GR . 

3. In all the above, the masses of light quarks and the effective mass of the gluon 
are neglected, whereas HERWIG uses these masses (only) in the final step, when 
it reconstructs the parton momenta from the values of x±,x 2 , cos 9. 

A. 1.3 Colour structure 

HERWIG treats all colour flows as distinct subprocesses. Where the colour flow is 
ambiguous, it is assigned according to the N — > oo limit. For example, in the FCR 
process gg — > QQ the colour-averaged matrix element squared is given by eq. (3.13). 
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In the N — > oo limit the two colour flows are proportional to u/t and t/u factors in 
the first bracket, and therefore we assign [16] 

da gg - da gg _ rt + l/ _ - x + _ 2/ _ 2 , da gg - x + _ 2/ _ 2 , (A.6) 

with <i<7 99 given in eq. (3.11). In dafg, which we shall call the t-flow contribution, 
gluon 1 is colour-connected to gluon 2 and the heavy quark Q, while gluon 2 is 
connected to 1 and Q. In d(jg a g \ the w-flow contribution, gluon 1 is colour-connected 
to 2 and Q, while 2 is connected to 1 and Q. 

In the other FCR process qq — > QQ, the colour connection is uniquely q — Q and 
q-Q. 



A. 2 Parton showers 

For a general introduction to the parton shower approximation, see Chapter 5 of 
ref. [17]. The parton shower in HERWIG [2] is performed using an angular variable 
£ [18] and energy fraction z. In the parton branching i — > jk, we have 

pj ■ p k Ej E k 

hijhik tLii tLii 

Thus £jk = 1 — cos 9 jk for massless partons. The evolution scale variable is Q = E^/^; 
thus the scale set by the above branching is Qi = Eiy/£j k . Colour coherence is 
simulated by angular ordering, which implies that the initial scales for showering on 
partons j and k are Qj = ZjQi and Qk = ZkQi respectively. 

The extension of the angular-ordered shower approximation to heavy quark pro- 
cesses is described in ref. [19]. 



A. 2.1 Initial conditions 

The initial conditions for parton showering on a given line % are determined by the 
invariant quantity Ef^ = p { - pj, where j is the colour partner of i. If i is a gluon line, 
it has two colour partners; one of these is selected at random, with equal probability. 
The showering on any line % is performed in a standard frame, in which the initial 
energy of i is E = Eij, its direction is along the z-axis while that of the colour 
partner j is in the (xz)-plane, and the upper limit on the angular evolution variable 
is £o = 1. Thus the initial scale for the shower is Q = E . 

In the FCR processes gg — > QQ, when the colour flow is the t-flow defined 
in app. A. 1.3 we have Eq = —t/2 for the outgoing heavy quarks, while for each 
incoming gluon we choose (separately) between Eq = —t/2 and Eq = s/2 with equal 
probability. When the colour flow is the w-flow, t is replaced by u. In the process 
qq — > QQ, we have Eq = —t/2 for all the parton showers. 
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A. 2. 2 Shower algorithm 

In all parton branchings, the leading-order massless splitting functions are used. 
In timelike (final-state) showering, the allowed region of z for the branching q — > 
qg is Q q /Q < z < 1 — Q g /Q where Q q = m^+VQCUT and Q g = m s +VGCUT, m q 
(^n,d,s[0.32,0.32,0.5]) and m 9 [0.75] being the quark and gluon effective masses, and 
VQCUT[0.48], VGCUT[0.10] minimum virtuality parameters. The argument of a s is 
z{\ — z)Q, which is an approximation to the k T of the emitted gluon 7 - see eq. (4.30). 
The corresponding Sudakov form factor is thus 

The next value Q' for the evolution variable is selected by solving A(Q') = A(Q)/R, 
where Q is the current value and R e [0, 1] is a random number. Branching stops 
whenever R < A(Q) is selected. The minimal value of Q is Q q + Q g . 

Spacelike (initial-state) showering is performed backwards, i.e. starting from the 
hard process. The evolution equation becomes 

A(Q')/f(x, Q') = [A(Q)/f(x, Q)]/R , (A.9) 

where / is the relevant PDF and x is the current value of the energy fraction of 
the spacelike parton. A compensating factor of f(x/z,Q) multiples the probability 
distribution in z. This "guides" the parton distribution to follow the input PDF. For 
example, since f(x/z, Q) = for z < x, branching to x values above 1 is prohibited. 

Perturbative branching of the spacelike parton in an initial-state shower stops 
when a value of Q < QSPAC[2.5] 8 is selected. However, if the parton is not a valence 
constituent of the corresponding beam particle, further non-perturbative branching 
is forced until a valence parton is generated. This is done in order to have a simple 
model of the beam remnant, composed of the other valence constituents carrying the 
remainder of the beam momentum. 

A. 2. 3 Azimuthal correlations 

After the shower, all the kinematics can be reconstructed from the values of z and £, 
except for the azimuthal angles. The distribution of these is isotropic unless dictated 
otherwise by the logical parameters AZS0FT[.TRUE.] and AZSPIN[.TRUE.]. If AZSOFT 
is . TRUE . , the azimuth of each branching is distributed according to the eikonal 
formula within the cone defined by the previous branching. In addition, if AZSPIN is 
also .TRUE., the azimuthal correlation in gluon branching due to gluon polarization 

7 The missing factor of y/2 is absorbed into the definition of the QCD scale A. 
8 The default value is relatively high because the input PDF paramctrizations may be unreliable 
(even negative) at lower scales. 



47 



is included. Polarization correlations for quarks are neglected, since they vanish due 
to helicity conservation in the massless limit. 

In the processes gi — > QQi, where % = q, q or g, there is an azimuthal correlation 
between the scattering plane and the plane of virtual gluon (g*) emission in initial- 
state branching, of the form [20] (in the collinear limit) 

M oc M gg P ( °\z) + C gg Q g *i(z) cos 20 , (A. 10) 

where is the azimuthal angle between the planes, is the leading-order % — > g 
splitting function, 

c -__ ^l^ + i-J-iV---'! (aid 

and 

Q gH (z) = -AC t , (A.12) 

with C g = C A = N and C q = C- q = C F = (N 2 - 1)/2N. 

Since this correlation vanishes in the limit m — > 0, it is also neglected in HERWIG. 
However, this means that there is matching between the parton showers and the 
matrix elements in the collinear limit only after azimuthal averaging of the latter. 
This poses a problem in MC@NLO implementation, whose solution is discussed in 
app. B. 

A. 2. 4 Momentum reshuffling 

After showering, each external parton line in the hard process has become a jet. 
The parton momenta have to be replaced by the jet momenta in such a way that 
energy-momentum is conserved, without seriously affecting the dynamics. This is 
achieved by the following momentum reshuffling procedure. 

Each final-state jet is first rotated from the standard showering frame (see 
sect. A. 2.1) to the direction of the corresponding parton in the hard process cm. 
frame, with a rotation about the jet axis to give the correct correlation with the di- 
rection of the colour partner. The magnitudes of the jet three-momenta in the hard 
process cm. frame are computed as follows. Suppose initially that each jet were given 
a three-momentum equal to that of the parton it replaces. Since the jet masses are 
not equal to the parton masses, this would violate energy-momentum conservation. 
However, energy-momentum conservation can be restored by rescaling the outgoing 
parton three-momenta in the hard process cm. frame by a common overall factor, 
to obtain the jet three-momenta in that frame. Once this factor has been computed, 
each jet is boosted along its axis to the required three-momentum. Note that the 
overall four-momentum of the hard process is not changed by this procedure. 

In the case of initial-state jets, the partons connected to the incoming beam 
particles are aligned with the beam directions (smeared by some intrinsic p T if that 
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is requested through the input parameter PTRMS[0]). Each of the two initial-state 
jets is then boosted longitudinally (i.e. along the beam directions) until the (now 
off-shell) momenta of the partons entering the hard process are consistent with the 
original kinematics. This is not a unique procedure, since not all kinematic variables 
can be restored to the values they had before showering. There are at present two 
options, controlled by the logical parameter PRESPL[.TRUE.]: 

• If PRESPL is .TRUE., then the cm. energy and longitudinal momentum of the 
hard process are preserved ("p-scheme"); 

• If PRESPL is .FALSE., then the cm. energy and rapidity of the hard process 
are preserved ("y- scheme"). 

Finally the transverse momenta of the initial-state partons entering the hard 
process are combined to give the transverse momentum of the hard process, and all 
the final-state jets are boosted transversely, by the amount required for transverse 
momentum conservation. 

In the previous HERWIG version (6.4) the longitudinal and transverse boosts 
were combined into a single boost in the direction of the new hard process momentum. 
However, this makes the matching of MC and NLO very complicated. Therefore the 
two boosts are performed separately in version 6.5. The difference between the 
two procedures corresponds to a rotation in the hard process cm. frame (Thomas 
precession). 

One can argue as follows that the effects of momentum reshuffling are beyond the 
next-to-leading order: to give at least one jet a mass requires one power of a s , and 
the average mass-squared is of order a s . Thus the fractional change in observables 
due to momentum reshuffling is expected to be of relative order a|. 

A. 3 Hadronization 

A. 3.1 Cluster formation 

After the perturbative parton showering and momentum reshuffling, all outgoing 
gluons are split non-perturbatively, into light quark- ant iquark or diquark-antidiquark 
pairs (the default option is to disallow diquark splitting). At this point, each jet 
consists of a set of outgoing quarks and antiquarks (also possibly some diquarks and 
antidiquarks) and, in the case of spacelike jets, a single incoming valence quark or 
antiquark. The latter is replaced by an outgoing spectator carrying the opposite 
colour and the residual flavour and momentum of the corresponding beam hadron. 

In the limit of a large number of colours, each final-state colour line can now be 
followed from a quark/ ant i- diquark to an antiquark/diquark with which it can form 
a colour-singlet cluster. By virtue of the preconfmement property of the shower [21], 
these clusters have a distribution of mass and spatial size that peaks at low values, 
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falls rapidly for large cluster masses and sizes, and is asymptotically independent of 
the hard subprocess type and scale. 

The clusters thus formed are fragmented into hadrons. If a cluster is too light 
to decay into two hadrons, it is taken to represent the lightest single hadron of its 
flavour. Its mass is shifted to the appropriate value by an exchange of 4-momentum 
with a neighbouring cluster in the jet. Similarly, any diquark-antidiquark clusters 
with masses below threshold for decay into a baryon-antibaryon pair are shifted to 
the threshold via a transfer of 4-momentum to a neighbouring cluster. 

Those clusters massive enough to decay into two hadrons, but below a fission 
threshold to be specified below, decay isotropically 9 into pairs of hadrons selected 
in the following way. A flavour / is chosen at random from among u, d, s, the six 
corresponding diquark flavour combinations, and c. For a cluster of flavour /1/2, 
this specifies the flavours f±f and ff 2 of the decay products, which are then selected 
at random from tables of hadrons of those flavours. The selected choice of decay 
products is accepted in proportion to the density of states (phase space times spin 
degeneracy) for that channel. Otherwise, / is rejected and the procedure is repeated. 

A fraction of clusters have masses too high for isotropic two-body decay to be a 
reasonable ansatz, even though the cluster mass spectrum falls rapidly (faster than 
any power) at high masses. These clusters are fragmented using an iterative fission 
model until the masses of the fission products fall below the fission threshold. In the 
fission model the produced flavour / is limited to u, d or s and the product clusters 
fif and // 2 move in the directions of the original constituents f\ and / 2 in their cm. 
frame. Thus the fission mechanism is not unlike string fragmentation [22]. 

In HERWIG there are three main fission parameters, CLMAX[3.35], CLP0W[2] and 
PSPLT[1]. The maximum cluster mass parameter CLMAX and CLPOW specify the fission 
threshold Mf according to the formula 

M™ = CLMAX™ + (mi + m 2 ) CLP0w , (A. 13) 

where mi and m 2 are the quark mass (RMASS) parameters for flavours f\ and / 2 . The 
parameter PSPLT specifies the mass spectrum of the produced clusters, which is taken 
to be M PSPLT within the allowed phase space. Provided the parameter CLMAX is not 
chosen too small, the gross features of events are insensitive to the details of the fission 
model, since only a small fraction of clusters undergo fission. However, the production 
rates of high-p T or heavy particles (especially baryons) are affected, because they are 
sensitive to the tail of the cluster mass distribution. Reducing CLPOW increases the 
yield of heavier clusters (and hence of baryons) for heavy quarks, without affecting 
light quarks much. For example, the default value gives no 6-baryons (for the default 
value of CLMAX) whereas CLP0W=1.0 makes the ratio of 6-baryons to 6-hadrons about 

1/4- 

9 Except for those containing a 'perturbative' quark when CLDIR=1 - see below. 



50 



There is also a switch CLDIR.fl] for cluster decays. For the default value, a 
cluster that contains a 'perturbative' quark, i.e. one coming from the perturbative 
stage of the event (the hard process or perturbative gluon splitting) 'remembers' its 
direction. Thus when the cluster decays, the hadron carrying its flavour continues 
in the same direction (in the cluster cm. frame) as the quark. This considerably 
hardens the spectrum of heavy hadrons, particularly of c- and 6-flavoured hadrons. 
CLDIR=0 turns off this option, treating clusters containing quarks of perturbative 
and non-perturbative origin equivalently. 

In the CLDIR=1 option, the parameter CLSMR[0] allows for a Gaussian smearing 
of the direction of the perturbative quark momentum. The smearing is actually 
exponential in (1 — cos#) with mean value CLSMR. Thus increasing CLSMR decorrelates 
the cluster decay from the initial quark direction. 

A. 3. 2 6-quark hadronization 

The process of 6-quark hadronization requires special treatment and the results ob- 
tained using HERWIG are still not fully satisfactory. Generally speaking, it is difficult 
to obtain a sufficiently hard B-hadron spectrum and the observed fr-meson/6-baryon 
ratio. These depend not only on the perturbative subprocess and parton shower but 
also on non-perturbative issues such as the fraction of 6-flavoured clusters that be- 
come a single B meson, the fractions that decay into a B meson and another meson, 
or into a 6-baryon and an antibaryon, and the fraction that are split into more clus- 
ters. Thus the properties of 6-jets depend on the parameters RMASS(5), CLMAX, CLPOW 
and PSPLT in a rather complicated way. In practice these parameters are tuned to 
global final-state properties and one needs extra parameters to describe fo-jets. 

A parameter B1LIM[0] has therefore been introduced to allow clusters somewhat 
above the Bit threshold mass M th to form a single B meson if 

M < M lim = (1 + BlLIM)M tft . (A.14) 

The probability of such single-meson clustering is assumed to decrease linearly for 
M th < M < M Um . This has the effect of hardening the B spectrum if B1LIM is 
increased from the default value. 

In addition, in HERWIG version 6, the parameters PSPLT, CLDIR and CLSMR 
have been converted into two-dimensional arrays, with the first element controlling 
clusters that do not contain a 6-quark and the second those that do. Thus tuning 
of 6-fragmentation can now be performed separately from other flavours, by set- 
ting CLDIR(2)=1 and varying PSPLT(2) and CLSMR(2). By reducing the value of 
PSPLT (2), further hardening of the B-hadron spectrum can be achieved. 

Figure 23 shows the b fragmentation function in Z° decay, i.e. the distribution 
of the energy fraction x E = 2E B /M Z in e + e~ collisions, where E B is the energy of a 
weakly-decaying 6-flavoured hadron in the Z° rest frame. The HERWIG predictions 
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Weakly decaying b hadrons 




Figure 23: Effect of HERWIG parameters on the b fragmentation function. 

for different values of some of the parameters discussed above are compared with 
the data of the SLD Collaboration [23]. For the MC@NLO predictions in sect. 8, 
parameters were left at their default values except for PSPLT(2)=0.5, corresponding 
to the long-dashed curve. However, as explained in sect. 8, all the results presented 
there are at the 6-quark level, and are therefore insensitive to the 6-hadronization 
parameters. 

B. MC subtraction terms 

In this section, we derive the MC subtraction terms <i£ a b| MC needed in order to define 
MC@NLO, eq. (2.1). We follow I, sect. A. 5. In the present case, the 0(af) term in 
the expansion of the result of HERWIG is given in eq. (5.1), which we report again 
here 



da 



MC 

ab L I 



= EEE rf <#" ■ (B-i) 

* ' * ' * ' MC 



This equation is identical to eq. (I. A. 58), except for the fact that more emitting legs 
and colour structures have been considered. 

We start from the case of initial-state radiation, considering emission from leg 1 
for definiteness. We can rewrite the relevant part of eq. (I. A. 58)) as follows: 



(+,0 



dxi i dx2ifi Hl \x 1 i/z^)f ( h H ' 2 \x 2 i)M l cb (x u ,X2i,t + ) d<j) 2 {xuX^S) 
S^i^'K^-d''). (B.2, 
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where we have explicitly indicated the integration over the azimuthal angle (p+ of the 
branching parton (which is trivial here but will not be so in the following). We have 
used the findings of sects. 4.1.2 and 4.2.2. The momentum fractions of the partons 
entering the 2—^2 hard process are x~u, and x 2i , given in eqs. (4.20) and (4.22) in 
the p- and y-scheme respectively (see sect. A. 2. 4). The t-channel 2^2 invariant is 
f+, given in eq. (4.17). The set (xu,x 2 i,t + ) fully specifies the momenta of the 2-^2 
hard process in the lab frame. In fact, the s-channel invariant is 

s = x u x 2i S = s+, (B.3) 

where s + is given in eq. (4.17), and the last equality holds thanks to eq. (4.18). With 
s+, t + , and u + = — s + — t + one gets immediately the 2 — > 2 momenta in the hard 
cm. frame: 10 

p 1)2 = £(l,0,0,±l) , h,2 = (E, ±kr, 0, ±h) , (B.4) 
where from eq. (3.7) 

E=\y/I^, k T = J _ m 2 ^ fe L = t+ , (B.5) 



Finally, the momenta in eq. (B.4) are boosted to the lab frame with 

yS = |iog?i. (B.6) 

In eq. (B.2), ^ s ^ ne Born matrix element squared, summed over spin and colour 
degrees of freedom (the average factor for initial-state ones is understood), times the 
flux factor, times a factor which depends on the partonic process and the colour flow. 
When cb = <p this factor is one; when cb = gg, it is one for / = s, whereas for I — t,u 
we have A4^ b = dagg / dfo, with da g l g given in eq. (A.6). Mil depends on xu and 
x~2i only through their product; however, we use the notation of eq. (B.2) order to 
remember the boost of eq. (B.6). 

We now perform some changes of integration variables in eq. (B.2). We replace 
(xu,X2i) with (x±,X2), the fractional parton momenta of the 2^3 process. We 
also replace (z^,^,ip + ) with (xi n ,yi n ,Lpi n ). Here, we do not need to define these 
variables precisely; they are used in order to eliminate the dependence upon the 
index / in the integration measure, and will be related to the three-body phase-space 
variables in the following. Finally, we use the fact that, in the cm. frame of the 
2^2 hard process, the two-body phase space can be written as follows 

#2(s) = ^dcos0 in , (B.7) 

107T 



10 In order to simplify the notation, the momenta have been rotated in the transverse plane along 
the direction corresponding to ip + = 0. 
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where the function f3(s) is given in eq. (3.8). We thus get 



d °ab 



(+,0 



^ab ,V) ®it' l) ^d>' l) dxi dx 2 dx in dy in dtp in d cos 6 in . (B.8) 



The factor C contains the dependence upon parton densities 

Cab ~ d{xi X2) jpfa (*«/*+)/«, (*»). (B.9) 

The actual form of the Jacobean that appears in this equation depends on the 
momentum reshuffling scheme (see app. A. 2. 4). By explicit computation, using 
eqs. (4.20) and (4.40) we find, in the p-scheme, 



d(Xu,X2i) xx. 



d(xi,x 2 ) y/x\ - x x x 2 v x v 2 ls 
In the ^/-scheme, we use eqs. (4.22) and (4.40) to obtain 

d{xiiiX 2 i) 



(B.10) 



d{x u x 2 ) 

The factor B collects the terms related to parton branching 

a 

e 



= x. (B.ll) 



-27r^ n ,y m ,^ n )27r *(0 U 1^ 



In contrast to the case of I (see eq. (I. A. 72)) the Jacobean has in general a non-trivial 
dependence upon the azimuthal angles; this is so in the cases I — t,u, and it is due to 
the complicated dependence of the HERWIG scale E upon the phase-space variables. 
Finally, the factor TC contains the full information on the hard 2^2 process, and 
phase-space factors 

n^ = ^^J^(xu,x x ,i + ). (B.13) 

We now turn to the case of final-state branching, considering emission from the 
heavy quark. The analogue of eq. (B.2) reads 



a<J ab 



= dx lf dx 2 ffi Hl \x lf )fl H2 \x 2f )M ab (x lf ,x 2f ,t Q )d(f> 2 (x lf x 2f S) 
x ^ d *W^^p(o>(4>)e ( ( J0)2 _ . (B.14) 

Here, the results of sects. 4.1.1 and 4.2.1 have to be used. The momentum fractions 
of the incoming partons are not affected by final-state branching, and eq. (4.7) holds; 
therefore 

s = x lf x 2f S = s = s Q . (B.15) 
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The t-channel invariant is given in eq. (4.5). The momenta of the hard 2^2 process 
are given in eq. (B.4) in the parton cm. frame, with E, k T , and k h computed as in 
eqs. (B.5) after the formal replacements (s + ,t + ,u + ) — * (sq, £q,mq)- These momenta 
are subsequently boosted to the lab frame by means of 

I/ffi = *log!^. (B.16) 

X 2 f 

Performing changes of variables analogue to those in eq. (B.2), we get the analogue 
of eq. (B.8) 



MC 



C W) B W) n W) dxi dx2 dXout dyout dVout d cos dout ? ( B . 17) 



where now 11 



4?'° = /i H,, 0Ei/)/r 2 W, (B.i8) 

BW , = i M^Sl ^«l e L$? - **) . ( b,9) 



ab ~ l&TT M abK x lf^ x 2f,tQ). (B.2U) 

Equation (B.8) and (B.17) are closely related to the MC subtraction terms 
dE ab \ MC /d(f) 3 needed for the definition of MC@NLO. We use eq. (I.A.61) 



d<T^% = dx 1 dx 2 ^- 

a<p3 



d<p 3 . (B.21) 



MC 



In order to proceed, we introduce an explicit parametrization for the three-body 
phase space. We use (see eq. (LA. 10)) 

dHs) = fjlf^ 1 ~ x)dx dy dtp d cos 9 Q , (B.22) 

where 8 Q is the scattering angle of Q in the QQ cm. frame. The variables x, y 
and (f refer to the emitted parton in the parton cm. frame of the 2^3 process: 
(1 — x) is the energy in units y/s/2, and y and tp are the cosine of the polar angle and 
the azimuthal angle respectively; x and y can be expressed in terms of invariants as 
shown in eq. (4.40). The ranges of the angular variables are < 9 Q < ir, < < 2ir. 
The parametrization of eq. (B.22) is identical to that of eq. (2.12) of ref. [3], provided 
that a (physically irrelevant) rotation of (p — 6 2 is performed there, in the transverse 
plane in the QQ cm. frame. The momenta parametrizations of eqs. (2.8) of that 



11 For final-state emission the Jacobean factor for the Bjorken x's is trivial, see eq. (4.7). 
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paper can thus be used here. We use eq. (B.22) to rewrite eqs. (B.8) and (B.17) as 
follows 



jv<+>') 



#3 

dKb 



4r B£» n[t' l) , (B.23) 
c (Q,i) B (Q,i) n (Q,i)^ (R24) 



MC 

where we have chosen (x in ,y in ,ip in ) and (x out ,y out ,(p out ) to coincide with (x,y,ip), 
and 12 

A+,l) _ d(xii,X2i) 1 f (H!)/- lAl)\AH2) ( - \ (T\o^\ 

L ab - *(„ - \ (I) J a \ X li/ Z +)h K X 2i), (B.25J 

o{x 1 ,x 2 ) z y 

Bca d(x,y,p) (R26) 

(+i0 _ 16tt d cos 9 in — { i) 

Kb ~ x lX2 Sdcos9 Q Mch{Xu ^ t+) ' (R27) 

4?'° = fi Hl) (*i)fl Ha) M, (B-28) 

*™ ^ f (^ _ ^) , (B , 9) 

Xi^-S 1 pfauXxS) dcos6 Q 

Here, we have used xs = XuX 2i S and s = x±fX 2 fS. The quantities xij and x 2i depend 
on xi, x 2 , and the invariants; thus, the Jacobean in eq. (B.25) must be computed at 
fixed (x,y,(p,9 Q ) = 3 . Furthermore, the Jacobeans in eq. (B.26) and (B.29) have 
to be computed at fixed 9 in and 6 out respectively. 

The MC subtraction terms appear twice in the definition of MC@NLO, eq. (2.1), 
since they contribute to the weights of EI and § events. These weights are treated 

(3) 

as ordinary MC weights in the MC evolutions, whose generating functionals are J-^c 
and J-^c respectively. Implicit in these generating functionals is also the dependence 
upon the initial conditions for the showers, i.e. the momenta of the 2^3 and 
2^2 processes. As the integral in eq. (2.1) requires, these momenta need to be 
specified for each point in the integration range (xi, x 2 , fa). This is straightforward 
in the case of H events, after a definite parametrization is chosen for the angles 9 in 
and 9 out , which allows the computation of the Jacobeans that appear in eqs. (B.27) 
and (B.30). 

The situation is more involved in the case of § events. Since here the initial 
condition for the shower is a 2 — > 2 process, a mapping of the three-body phase space 

12 Here and in what follows we have some abuse of notation, since we denote by the same sym- 
bols C, B, and H quantities which possibly differ by multiplicative factors from those previously 
introduced. 
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onto the two-body one is necessary. This map is provided by the MC, which relates 
the momenta of the partons entering the 2 — > 2 hard process to those of the partons 
emerging from the shower. In the case of a single branching, this map has been 
denoted by Vu^n in sect. 2. In this section, we have given explicit implementations 
of Vu^s- in the case of initial-state emission, one uses eqs. (B.4)-(B.6), where the 
2^2 invariants are expressed in terms of the 2^3 ones as explained in sect. 4.1.2. 
For final-state emission, analogous relations hold (see the discussion after eq. (B.15)). 
These maps also provide us with the definitions of the scattering angles 0i n and 9 out 
in terms of the three-body phase-space variables. In fact, the t-channel invariant is 
always given by the function (see eq. (3.9)) 



s and 6 being the s-channel invariant and scattering angle in the hard cm. frame 
respectively. Thus, taking into account eq. (4.8), and the fact that sq = s, we have 



which can be solved for 8 in and 6 out using eqs. (4.17) and (4.5) respectively. 

However, eqs. (B.32) and (B.33) manifestly define two different maps and 



whereas eq. (2.1) requires p£3 = = Vm^s- In fact, if p£S ^ V i 



for a given three-body configuration (xi,x 2 ,03) one would get two two-body config- 
urations generated by V^l s and V^l, and the choice of which one to use in as 
initial condition for the shower would result in an ambiguity in MC@NLO. Besides, 
the complicated functional relations implicit in eq. (B.32) and (B.33) would make the 
task of the implementation of event projection in the NLO cross section a difficult 
one. 

A solution for the first problem mentioned above is to treat initial- and final- 
state emissions independently in MC@NLO. The NLO cross section must also be 
written accordingly, since MC subtraction terms act as local counterterms. This can 
be achieved for example using the technique of ref. [20], which basically amounts to 
partitioning the phase space into regions dominated by initial- or final-state emis- 
sions. 13 In terms of numerical accuracy and unweighting efficiency this may be the 
best strategy in those cases in which final-state collinear emissions are singular at 
the level of short-distance cross sections, such as in jet physics. However, it would 
still leave us with the problem of an easy implementation of event projection. 

We argue that, regardless of the presence of final-state singularities in real matrix 
elements, the simplest possible definitions of the maps and should be 

adopted. We now show how this can be achieved; as a by product, we also show 

13 Ref. [20] is based on the subtraction method, so no approximation is involved in this phase-space 
partition. 



t(s,9) 



±s(l-/3(s) cos#), 




t(s + V!+ v 2 ,9 in ) = t 
t(s, 9 out ) = tQ , 



(B.32) 
(B.33) 
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that, in the present case, initial- and final-state emissions can indeed be treated 
simultaneously. 

We start by observing that momentum reshuffling is an effect beyond NLO (see 
sect. A. 2. 4). It follows that two 2 —» 3 kinematic configurations are equivalent at the 
NLO if they coincide in the limit in which the off-shell parton (i.e., the parton which 
has branched) has its on-shellness restored. This is a consequence of the fact that the 
branching procedure is exact only for zero-angle emission. Thus, for a given 2^3 
configuration we can take its zero- angle-emission limit, and extract from there the 
2^2 momenta we seek. This is most easily done in the partonic cm. frame, where 
eq. (B.4) holds; eqs. (B.5) are then used to obtain E, k T , and k L , with the 2^2 
invariants expressed in terms of the zero-angle limits of the 2^3 invariants. The 
limits can be computed by choosing any parametrization for the three-body phase 
space; in the following, we shall use eq. (B.22), exploiting the explicit formulae of 
ref. [3]. 

For initial-state emission (say, from parton 1), from eq. (4.17) we get, in the 
zero-angle-emission limit 

-(o) 

s + = xs , 

^ = -\ xs [l - p(x S ) cos# Q ] . (B.34) 
Replacing t+ with in eq. (B.32), and solving for 6 in , we get 

cos 6 in = cos 0q ■ (B.35) 

This is what we expect, since for strictly collinear initial-state emission, the parton 
2^2 cm. frame (where 6 in is defined) and the QQ cm. frame (where 9 Q is defined) 
coincide. 

For final-state emission from the quark leg, from eq. (4.5) we get 

7(0) 



where 



Note that 



*Q = 'I 3 i 1 - M cosd out(x, cos9 Q )] , (B.36) 



„ / „ x 1 — a; — (1 + x) cos^o „„ s 

cos9 out (x,cos6 Q ) = -— }- B.37 

1 + x — (1 — X) COS# Q 



cos 0^(1, cos Q ) = cos6» Q , (B.38) 

which is what one expects on physical grounds, since x — ■> 1 is the soft limit in the 
parametrization of eq. (B.22), and it is only in this limit that the parton 2^2 cm. 
frame and the QQ cm. frame coincide in the case of final-state emission. 

We can now use eqs. (B.35) and (B.37) in eqs. (B.27) and (B.30) respectively. 
This fully defines the MC subtraction terms, at least in the case of H events, where 
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the initial condition for the MC shower is given by the 2^3 configuration, as 
specified by the integration variables (x\,x 2 , 03). However, we still have the problem 
that 7^ T^w^l- To see this explicitly, we write the 2 — > 2 momenta that we 

obtain with the zero-angle-emission prescription: 

• for initial-state emission 

Pi,2 = (1,0,0, ±1) , h,2 = ±Vx^(l,±P(xs)sm6 Q ,0,±p(xs)sm6 Q ) ; (B.39) 

• for final-state emission 

Pi,2 = |«s/i(l,0,0,±l) , h,2 = ^(l,±P(s)sine out ,0,±P(s)sme out ) , (B.40) 

with 9 out given in eq. (B.37). The momenta in eq. (B.39) are then boosted to the 
lab frame using eq. (B.6), those in eq. (B.40) using eq. (B.16). We now observe that 
the momenta in eq. (B.39) can be obtained with the parametrization of eq. (B.22) 
by imposing the soft limit x — > 1, but freezing the cm. energy to the value xs. 
Furthermore, as described in sect. 4.1.1, for final-state emission the reshuffling does 
not change the direction of the outgoing heavy quark in the partonic cm. frame; 
thanks to eq. (B.38), a practical way of computing it is to consider the soft-emission 
limit in the parametrization of eq. (B.22). Thus, in the soft limit the scattering 
angles in the parton 2 — > 2 cm. frames coincide for initial- and final-state emission. 
Unfortunately, the 2 — > 2 momenta do not coincide, because of the different s-channel 
invariants and boosts to the lab frame. However, we can force them to coincide with 
a simple formal manipulation, similar to what is done for event projection. We can 
use the identity xs = XuX 2 iS in eq. (B.39), and s = X\X 2 S in eq. (B.40). Next, we 
formally replace x± with x u and x 2 with x 2i in e q. (B.24) (which is allowed, these 
variables being just integration variables there), and then we change integration 
variables back to x\ and x 2 . Eqs. (B.28) and (B.30) become 

cT = ^il^/W^f*)^), (B.41) 

0{Xi, x 2 ) 

0,(0,1) 167r dcos9 out —(i) _ 

= x u x 2t S dcos6 Q M <* { ^ X ^ tQ) ■ (R42) 

In may appear counterintuitive to have a dependence on x H and x 2i in these equa- 
tions. However, it must be clear that this has nothing to do with initial-state emis- 
sion. Simply, the freedom of changing integration variables allowed us to scale the 
partonic cm. energy in eq. (B.24) by a factor of x. Thus, eq. (B.39) now holds 
for final-state emission as well. Furthermore, the boost to the lab frame is now 
performed with eq. (B.6) rather than with eq. (B.16). This is precisely what we 
wanted to achieve, since a unique Vu^s has now been defined. For a given choice of 
(xi,X2,</>3), one reconstructs the 2^2 momenta in the hard cm. frame by taking 
the soft limit with the s-channel invariant fixed to xxiX 2 S, and then boosts to the 
lab frame using eq. (B.6). 
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We stress that, in the case of final-state emissions, we define the MC subtraction 
terms using eqs. (B.28)-(B.30) for EI events, and eqs. (B.41)-(B.42) for § events. 
In principle, we should therefore introduce the notations dE^\ MC and cE^| MC , but 
we prefer to avoid it, since it has been shown that the difference between the two 
amounts to a change in the integration variables. 

The zero-angle procedure implies the possibility of several definitions of the MC 
subtraction terms, which differ beyond NLO. For example, one could use the 2-^2 
invariants obtained in the zero-angle-emission limit, eqs. (B.34) and (B.36), in the 
computation of the hard 2^2 matrix elements, for both EI and § event, or only for § 
events. We prefer not to consider the latter option, which results in a total rate com- 
puted by the MC@NLO different from the NLO one. The former choice has indeed 
been implemented; the results have been found to be identical to those obtained with 
our default choice in the case of top production. When bottom production is con- 
sidered, very small differences are visible in the tails of those distributions which are 
effectively of leading-order accuracy in the fixed-order computation of ref. [3], such 
as the transverse momentum of the pair or A(f>. However, these differences never 
exceed a few percent, and are much smaller than the uncertainties of the fixed-order 
results due to scale variation. This confirms that reshuffling effects are beyond the 
accuracy of MC@NLO. 

As pointed out in I, sect. A. 5, the MC subtraction terms given in eqs. (B.23) 
and (B.24) cannot act as local counterterms for real emission matrix elements, since 
the angular distribution of a soft gluon emitted by the MC does not agree with 
the corresponding perturbative result. However, it has been argued in I that this 
effect must be irrelevant in the definition of observables, at least for infrared-safe 
ones. Thus, in I the MC subtraction terms were defined by smoothly matching 
what is obtained from the perturbative expansion of the MC result with the leading 
singular behaviour of the real matrix elements in the soft limit. The matching has 
been defined through a parameter-dependent damping function, eq. (I. A. 86). The 
physical observables were found to be independent of the parameters used to define 
the damping function, thus practically confirming that infrared-safe observables are 
insensitive to the angular distribution of soft-gluon emission. 

Here, we shall follow the same strategy. Since the soft singularity structure is 
more complicated than in the case of gauge-boson pair production, and in particular 
cannot be associated with the soft divergence of an Altarelli-Parisi splitting function, 
we generalize eq. (I. A. 83) and eq. (I. A. 84) as follows: 



ah 



(L,l) 



= Gs(x)g c (y) J2 E ~it~ + 0- - Os(x))M ab (S) 
i —~ a(p3 mc 

+ (1 - g c (y))M ab (C) - (1 - g s (x))(l - Q c (y))M ab (SC), (B.43) 



where M a b(S), M-ab{C) and M a b(SC) denote the leading singular behaviour of the 
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real matrix element in the soft, collinear, and soft-collinear limits respectively. The 
functional form of the damping functions Q s (x) and Q c {y) is the same, and is given 
in eq. (I. A. 86) (with x DZ = 0); the subscript is to remind us that the parameters 
these functions contain can be varied independently. The function Q s (x) smoothly 
approaches zero in the soft limit 

lim = , x = 1 + - + V2 , (B.44) 

x-*l S 

and the function Q c (y) smoothly approaches zero in the initial-state collinear limits 

lim C c (y)=0, y = ^^l. (B.45) 

Note that the function Q c (y) would not be necessary away from the soft limit, if 
azimuthal correlations (see app. A. 2. 3) were either absent or properly described by 
the MC. As in I, no evidence has been found of a dependence of the physical results 
upon the parameters entering the damping functions. 

In summary, eq. (B.43) is our final formula for MC subtraction terms. The 
quantities cE^'^| MC are defined in eqs. (B.23) and (B.24) for initial- and final-state 
emissions respectively. C, B, and 7i appearing there are given in eqs. (B.25)-(B.27) 
for initial-state emissions. For final-state emissions, we use eqs. (B.28)-(B.30) for HI 
events, and eqs. (B.29), (B.41), (B.42) for § events. 



C. Colour flow codes 

We list in tables 3 and 4 the codes IC used in MC@NLO to transmit colour flow in- 
formation to the HERWIG event generator. The convention is that q and q represent 
the colour and anticolour of parton i. If i is a quark (antiquark) then its anticolour 
(colour) is zero. Thus for example when IC = 1, corresponding to qq — > gQQ (see 
table 4), we have c q = c g = 1, Cq = cq = 2, c g = cq = 3, meaning that q and 
g are colour-connected, q and Q are anticolour-connected, and the anticolour of g 
is connected to the colour of Q. In accordance with the Les Houches convention, 
the non-zero colour labels entered into the event common block HEPEUP are actually 
IC0LUP(1, i) = 500 + a and IC0LUP(2, i) = 500 + q. 

D. £ subtraction 

In I, an NLO subtraction scheme was introduced, called ( subtraction, which al- 
lows one to reduce the number of negative-weight events occurring in MC@NLO, 
compared to that resulting from the implementation of the "standard" subtraction 
formulae of ref. [24], in which the computation of W + W~ cross sections to NLO 
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IC 


12 - 


-> 34 


CiCi 


c 2 c 2 


c 3 c 3 


C4C4 


1 


qq- 


*QQ 


10 


02 


10 


02 


2 


qq - 


*QQ 


01 


20 


20 


01 


3 


99 ~ 


+ QQ 


12 


23 


10 


03 


4 


99 - 


+ QQ 


12 


31 


30 


02 



Table 3: Colour flow codes for 2 — > 2 configurations 



IC 


12 - 


-> 345 


C1C1 


C2C2 


C3C3 


C4C4 


C5C5 


1 


QQ — 


+ qQQ 


10 


02 


13 


30 


02 


2 


QQ — 




10 


02 


32 


10 


03 


3 


QQ — 




10 


21 


30 


20 


03 


4 


QQ — 




10 


23 


20 


10 


03 


5 


QQ — 


* #g<2 


01 


20 


23 


30 


01 


6 


qq - 


+ </QQ 


01 


20 


31 


20 


03 


7 


qg - 


^<?gg 


01 


23 


03 


20 


01 


8 


Q9 ~ 




01 


12 


03 


30 


02 


9 


gq. - 




12 


20 


30 


10 


03 


10 


gq - 




12 


30 


10 


30 


02 


11 


gq- 


^<?gg 


12 


03 


02 


10 


03 


12 


gq- 


+ <?gg 


12 


01 


03 


30 


02 


13 


gg - 


^gg 


12 


23 


14 


40 


03 


14 


gg - 


+ <?gg 


12 


34 


32 


10 


04 


15 


gg - 


^gg 


12 


23 


43 


10 


04 


16 


gg - 


^gg 


12 


31 


34 


40 


02 


17 


gg - 


^gg 


12 


34 


14 


30 


02 


18 


gg - 


-gQQ 


12 


31 


42 


30 


04 



Table 4: Colour flow codes for 2 — > 3 configurations 



accuracy was originally performed. The ( subtraction is essentially denned by the 
requirement that the subtraction terms be non-zero only in the region 

^<C, (D.l) 
s 

k T being the transverse momentum of the real parton emitted at the NLO. In terms of 
the variables x and y (see eqs. (B.44) and (B.45) for their boost-invariant definitions, 
and eq. (B.22) for an explicit phase-space parametrization which uses them), eq. (D.l) 
becomes 

V(x,y) ee (l- X y(l-y*) < C. (D.2) 
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As pointed out in I, eq. (D.2) is not specific to vector boson pair production, and 
can be applied to any NLO cross section. However, the formulae for the partonic 
cross sections derived in app. A. 3 of I need to be generalized. That is the aim of this 
section. Although we shall closely follow the notation of ref . [3] , it will be clear that 
our results are relevant to any hard production process of a heavy (or hard) system, 
for which the matrix elements have all possible soft and collinear divergencies, except 
those due to final-state collinear emissions. 

We start by observing that eq. (D.2) implies that there is a minimum value of 
x, for fixed y, equal to 



p = mm(x) = max [ 1 - ^/ T~~2 > P ) > ( D - 3 ) 

where p is the minimum value of x allowed by kinematics (for heavy flavour produc- 
tion, p = 4m 2 / s). We then write the analogue of eq. (3.5) of ref. [3] in a more general 
form 14 

da$ = d$W^dxdyd<p (I - x)- 1 - 2 * (I - y 2 )- 1 -* \ m ab (x,y,v) , (D.4) 

where d<&£^ is the (i-dimensional (d = 4 — 2e) reduced phase space for the production 
of the heavy system with a squared invariant mass equal to sx. In d$r^ we also 
incorporate further normalization factors, that reduce to 1 in 4 dimensions. Thus, 
the phase space of the heavy system plus the light parton is 

d$ = d<$>W-—dxdydip (1 -xf- 2t {l -y 2 )~ e \ sin <p\~ 2e . (D.5) 

647T d 

(x) 

In the present context, the precise (i-dimensional form of d& r is irrelevant. The 
only thing we need to know is that in 4 dimensions it is equal to the phase space of 
the reduced system. Thus, in the case of heavy flavour production 

^U = 4 = %^cos0 Q , (D.6) 

107T 

which coincides with eq. (B.7) for s — > xs. Using this equation in eq. (D.5), we get 
d<3> = d(f) 3 (s) in 4 dimensions, with d(f) 3 (s) given in eq. (B.22). In eq. (D.4), we have 
also defined m a b as 

m ab = (1 - x) 2 (1 - y 2 ) M ab (x, y, (p) , (D.7) 

where A4 a b(x,y,(p) is the invariant spin-averaged squared amplitude divided by the 
flux factor. Its dependence upon the kinematic variables of the reduced system (9 Q 
in the case of heavy flavour production) is not explicitly shown. 



14 



Consistently with app. B, and unlike the convention of ref. [3], we have here < ip < 2ir. 
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Following the reasoning of ref. [3], we now use the expansion 

'log(l -x) 



1 — x 



-l-2e 



B~ Ae ( 1 

-*(!-*)+( — 



2e 



x 



It 



1-x 



+ 0(e 2 ), (D.f 



with (3 = y/1 — p (notice that here p depends on y, which is not the case in ref. 
This yields 

,1-e 



d°$ = + dQM^dxdydip (1 - y 2 )' 1 ^ | sin^| 



-2e 



X 



1 — X 



2e 



log(l — x) 
1 — X 



m ab (x,y,ip), 



(D.9) 



where 



da2 = d$ r ^dyd<p(l-y 2 ^ - 



sin ip 



-2c 



64tt 3 

with d$ r = d$^\ x= i. We can now expand 



(3- 



-4c 



2e 



m ab (x,y, <p)\ x= i , (D.10) 



(i - y 2 )- 1- * = + 1/) + *(i - y)] 

where we define 



(2cD 



1 

2e + 2 



+ 



1 -yJs> \ 1 + y 



u = 1 — \ max ( 1 — 



c 



(1-x)' 







+ 0(c), 
(D.ll) 

(D.12) 



This quantity has the same formal meaning as u in ref. [3]; unlike iv, however, it does 
depend on x. Using eq. (D.ll) in eq. (D.9) we get 



da (r) = j («) , j (c+) , . (c-) , , (/) 

a0 ab UO ab ^ UO ab ^ U0 ab ^ U0 ab ' 



(D.13) 



where 



A ( c± ) 
d <b = 



,1-e 



— — dx dip I sin 09 1 2e 
r 64tt 3 1 1 



1 — X 



2e 



log(l — x) 
1 — X 



X 



2e 



+ 



1 logo; 
21 -x 



m ab (x,y,ip)\ 



y=±l , 



and 



da 



(/) 

ab 



2 Vl-x 



m ab (x,y,ip) 
1-x 



(D.14) 



(D.15) 



Notice that p prescriptions, rather than p prescriptions, appear in eq. (D.14), since 
p = p for y — > ±1. Furthermore, the expression logu)/(l — x) in eq. (D.14) does not 
need a regularization prescription, since as x — > 1 also cD — > 1. 
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( f) 

To see precisely the structure of the generated counterterms in da^ b , let us 
expand the expression 

dxd y(-rh) ~(rb) v{x ' y) (D ' 16) 

according to the definition of the distributions. Here V(x, y) represents symbolically 
some function of the final-state kinematics, regular for x — > 1 and y — > 1. We get 

(D.17) 



1 -^/s V 1 -y 



dxdy r y) - e(o) - (1 - 1) G(x - p)V(l,y) - e(y)V(l, 1) 



1 — x I 1 — y 1 — y 

where in the middle (last) expression we have used the fact that u — > 1 when x — > 1 
(p — > p when y — > 1). Using the relation 

e(u-(l-y)) = e(x-p)G(y) (D.18) 

we can rewrite the above expression as 

1 



J dxdy 



x J p\ l ~y 



V(x,y) 



dxdy / V(x, y) e(x - p)(V(l, y) + 0(y)(V(x, 1) - V(l, 1))) 



/" dxdy j V(x,y) _ 
J l-x\ l-y 



+ e(p-,)e ( y)y(M) 



i-y 

The last term is finite by itself, since 

&(p-x) = &((l-x) 2 (l-y 2 )-(). (D.20) 
It is a soft term, which we can rewrite by explicitly computing the integral over y 

/•^e(p-*)e W y(i,i) /•■^ (D21) 

J 1 - x l-y J p 1 - x 

For consistency with I, we include this term in da^ b + \ which is then modified with 
the replacement 

lDg " - l0g "-5(l-x) fte*** = . (D.22) 



1 — x 1 — x J p 1 — x \ 1 — x / fi 
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We now present the corrections to be added to the soft and collinear terms in order 
to go from the standard subtraction scheme of ref. [3] (computed with min(x) = p 
and ou — 1) to the ^-subtraction scheme. In the soft term, the correction is given by 



,1-e 



e 

x / dydtpil-yY^lsmtpl-^ip-^ ~ (3-* e )rn ab (x,y,<f)\ x=1 . (D.23) 



-2e th-U 



We find 



= e(i- y 



c 



i-y 2 

2 C" 



c 



P 4 



i-y 2 

C//9 4 



-4c 



(D.24) 



Since the function implies \y\ < 1, we only need to expand the square bracket to 
order e, thus getting 



/r 4e - p- 



Ac 



e[i-^-r 



log — + log(l - y 2 ) 



(D.25) 



So, our final formula for the correction is 



(«) 

ah 



128tt s 



X 



/ dip 



log^ | log(l-y 2 

X _ y2 1 — ?/ 2 



m a6 (a;,y,</?)| 



x=l 



with (see eq. (I.A.27)) 



i/ = e(i-i 



(D.26) 



(D.27) 



The soft limit m a b(x, y, (p)\ x =i is non-zero only if the radiated parton is a gluon. It 
is given in general by a sum over eikonal factors 



m ab (x,y,ip)\ x= i = y2(qi,q m ) C Xm , (qi,q m ) = — 

' * Hi • 



qi ■ q m 



Im 



qi-kq m -k' 



(D.28) 



where q\i m are the external momenta of the reduced process, k is the momentum 
of the soft gluon, and Ci m are functions of the momenta of the reduced process. 
Expressions for m ab (x, y, (p)\ x =i in the case of heavy quark production can be easily 
obtained (after correcting a couple of misprints 15 ) from Appendix A of ref. [3]. Terms 

15 1/2C A in the last line of eq. (A. 12) should read C A , and the —2 in the last line of cq. (A. 21) 
should read +2 
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arising from the eikonal factor (pi,p2) associated with initial-state soft emissions can 
be integrated analytically over y and tp in eq. (D.26). For the other terms, only 
the azimuthal integration can be easily computed analytically. The y integration is 
performed numerically. 
The collinear correction is 

Ada^^d^^dxdip^^j m ab (x,y,<p)\ y=±1 . (D.29) 

Finally, for the finite part of the real cross section an equation identical to eq. (I. A. 26) 
holds: 16 



+ 



[l- X ){l-y)) T y(l- x )(l + y )) v 

where the P-distribution prescription is defined by 



m *( x >V>ri d * , (D.30) 
1 — x 



, , . , , V(x,y)= (D.31) 

(1 -x)(l ±y)J v 

1 f V(x, y) Q(x - p)(V(l, y) - Q(Ty)(V(x, fl) - V(l, Tl))) ] 
l-x \ l±y l±y J ' 

We now discuss the relation of the above results with the results of I. The 
collinear correction is equivalent to the term 

log(l- .T c (aO) 

X / p 

of eq. (I. A. 25). In the case of the soft correction, only the terms corresponding 
to {l,m} = {1,2} are non-zero in eq. (D.28); these can be integrated analytically, 
yielding the T s term in eq. (I. A. 24). 

In the case of heavy flavour production, we find that ( subtraction reduces only 
marginally (l%-2%) the number of negative weights with respect to standard sub- 
traction. In fact, due to the simultaneous presence of several colour flows which 
induce different dead zones (see figs. 2 and 3), the subtraction region of eq. (D.2) 
never matches closely the region in which HERWIG radiation is allowed. Nevertheless, 
it is still advantageous to use ( subtraction, since the parameter tuning, necessary 
to reduce as much as possible the number of negative weights, is easier than in the 
case of standard subtraction (simply because the latter depends on two parameters 
rather than one). 



16 In eq. (I. A. 26) a different notation is used, which leads to an error if one forgets to freeze P(x, y) 
when applying the V prescriptions, as required in eqs. (I.A.21)-(I.A.23). 
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